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CHAPTER I 

INTRODUCTION 


Satisfactory performance of the machines can be 
ensured only through rigorous and accurate analysis and 
synthesis of mechanisms involved in the system. Improper 
design of mechanisms may lead to loss of accuracy and also 
eventual failure of its members. Under normal circumstances, 
the analysis and synthesis of the mechanisms are done without 
considering the finite rigidity of the members and the 
problem is mainly geometric. However, in many cases the 
design of the members also involves the consideration of the 
static and dynamic forces acting on the mechanism. ¥/hen 
the speed of operation is low, the dynamic forces play a 
relatively insignificant role. But as the speed of the 
mechanism increases, the accelerations of the different 
members of the mechanism and consequently the dynamic 
forces increase. For sufficiently large speeds, the dynamic 
forces may play even the predominant role. Under such 
circumstances, the members ai^e designed from the consideration 
of the stresses developed due to the dynamic and static 
forces. The usual practice, until recently, had been to 
consider the links to be absolutely rigid. This is 
acceptable so long as the accuracy required in the 
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relationship between the input and the output is not very 
high. Apart from this, in high speed mechanisms, the 
chances of the mechanism being subjected to resonant condition 
is quite high. In the last decade, machines comprising of 
high speed mechanisms have been developed to meet the demand 
of the industry. Moreover, the accuracy requirement of such 
machines is also having an increasing trend. The elastic 
displacements, which are the consequences of increased 
speeds and loads, may cause inaccuracies of position in 
addition to noise and fatigue. For such cases, the usual 
design procedures for the mechanism may not yield satisfactory 
results . 

In view of the above facts, the research work in 
the field of analysis and design of the mechanisms with 
flexible members is gaining momentum. Already a considerable 
amount of work in this area has been published. ITevertheless , 
the extent of research in this field is far from complete 
and enough scope is present to solve a number of challenging 
problems. 

Introduction of the flexibility of the links 
complicates the analysis and synthesis of the mechanisms 
considerably. Analytical treatment of such problems leads 
to a number of extremely complicated differential equations 
solution of which, even by numerical methods, requires very 



simplifying assximptions . To overcome such difficulties, a 

number of numerical methods have been proposed by different 

workers. However, such methods have neglected many important 

dynamical effects whose contributions may be significant. 

Therefore it has become necessary to assess the magnitudes 

of these effects. Moreover, efforts should be made to evolve 

new techniques to solve such problems with reduced computer 

time without detoriating the accuracy level. Apart from 

these, a mechanism with elastic members should be also 

studied to investigate the stability of the system. 

Before starting writh the main body of the thesis, 

it is pertinent to review the extent of research published 

to date in the fields of the flexible mechanisms studied in 

the present work. A survey of these investigations (and 

many others also) is carried out in the two excellent review 

( 1 2 ) 

papers ’ published a few years ago. In the next section 
only the important ones with recent origin are reviewed 
briefly. 

1.1 Review of Past Works 

1.1.1 Deflection analysis 

( 5 ) 

In 1965j Kozsevnyikov^ reported on the flexural 
vibrations of an elastic coupler of a slider-crank mechanism 
considering only the effect of the transverse accelerations. 


i 


Ueubauer, Cohen and Hall^'^\ in the analysis of the same 
mechanism, considered the effect of a constant axial force 

(s') 

derived from the rigid body analysis. Viscomi and ' 

contributed several improvements on the analysis of the same 
mechanism. The authors include the Coriolis components of 
acceleration as well as the additional relative normal and 
tangential components of accelerations arising from the 
transverse vibrations of the coupler in an acceleration 
field. The effect of the axial forces due to the rigid-body 
normal acceleration as well as due to the foreshortening of 
the coupler is also taken into account. The further simpli- 
fication of the differential equations excludes all effects 
of longitudinal vibrations and retains only the first two 
modes of the transverse vibrations. Jainiski, Lee and Sandor 
found out both the steady-state axial and transverse 
deflections in the analysis of the same mechanism without 
the piston force and with the piston force But 

subsequent linearisation of the differential equations 
neglects the effect of all axial forces on the bending moment 
of the coupler. Furthermore, the solution method requires 
that the obliquity ratio must be small enough to allow an 
asymptotic expansion of the deflections in terms of this 
ratio. This restricts the validity of the analysis to 
limited cases only. 



( 8 ) 

Meyer zur Gapellen^ ^ analysed the transverse 
vibrations of an elastic coupler of a four-bar mechanism 
which included the effect of the rigid body axial forces 
only. But the analysis is made simple by equating the 
periodic terms of the axial forces to zero and the final 
solution is rather approximate. 

(q) 

Using a lumped mass model, Sadler and Sandor ^ 
analysed the slider-crank mechanism and also for the case 
when the coupler has a concentrated mass at the end of its 
rigid extension. The method of analysis is later extended 
to a planar four-bar linkage composed of elastic links 
The longitudinal deformations resulting from the axial forces 
and the foreshortening of the links due to bending are 
neglected. The investigation is confined to the transverse 
deflections only by employing D' Alembert's principle for 
each lumped mass and solving the Euler's equation of the 
elastic curves by finite difference method. The additional 
relative normal and tangential accelerations are included 
and the pin forces are calculated from the force and moment 
balance of each link at their instantaneous positions. In 
addition to the gross rigid body motions, the analysis 
takes into account the additional rigid body motions of the 
coupler and the follower due to elastic deflections of the 
crank. Because of the assumption of the inextensibility of 
the links, the effects of the elastic deformations of each 
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link on the rigid body motions of the other can not be fully 
studied by their analysis. It is demonstrated through an 
example that the rigid body motions of the follower due to 
the elastic motion of the crank is comparable to its elastic 

motions and should not be neglected. 

(11 12 ) 

Winfrey ’ first used a finite element model 
to analyse the longitudinal, transverse and torsional dis- 
placements of , a general mechanism. Six elastic and three 
rigid body co-ordinates are defined for each link (planar). 
This type of modelling of the links increases the size and 
complexity of the matrices without any corresponding increase 
in the accuracy of the displacements. Kinematic analysis is 
performed by adopting the general iterative procedures of 
Denavit and his group ^ . A connectivity matrix is 

constructed to assemble the element matrices into the system 
matrices. Since the connectivity matrix is to be constructe|i 
for each problem separately, it is seldom preferred to 
other techniques of assembly available in the finite elemen 

i 

method. The author directly applies the principle of the 
conservation of generalised momentum during the free vibra 
of the mechanism to form a constraint equation for the 
elimination of the rigid body degrees of freedom. Here 
again the connectivity matrix and its derivatives are 
involved in addition to numerous calculations. Upon 
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utilisation of this constraint equation in a transformation 
to the generalised co-ordinates, the stiffness matrix is 
made non-^singular and the eigenvalue problem in the trans- 
formed co-ordinates is solved by inverting this non-singular 
matrix. The final transformation of the original equation of 
motion to the normal co-ordinate system is achieved by using 
first the constraint equation and next the modal matrix 
formed after solving the eigenvalue problem. Since the 
constraint equation is valid for the free vibration only, 
the extension of its range of applicability to the equation 
of forced vibration is unwarranted. The equations in the 
normal co-ordinate system are made, uncoupled by assuming 
the damping as proportional to the mass matrix and solved by 
Laplace transform method. The final values of the displace- 
ments and the velocities found out in the system co-ordinates 
of the current interval of time are used as the * initial 
conditions' for the system co-ordinates of the next interval. 
This is justified as the system co-ordinates are located 
at element positions for each interval. 

The formation of the system force vector uses the 
connectivity matrix and its derivatives alongwith the unknown 
displacements. Thus the equations of the forced vibrations 
are to be solved by iterative procedures. In an initial 
analysis, the author finds that the inertia forces due to 
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the relative accelerations are only 1% of the inertia forces 
arising from the rigid body accelerations and these terms 
are dropped out in the subsequent calculations as their 
retaintion entails intricate calculations within the iteration 
loop. This finding may not be true in general. Under another 
assumption, the force vector term is reduced to its final 
form which does not take into account the variation of the 
accelerations along the links. The effect of the axial 
forces on the transverse vibrations of the links is altogether 
neglected in the above analysis. 

(15) 

Erdman, Sandor and Oakberg analysed mechanisms 
by the flexibility approach of the finite element method. 
Starting from the rigid body inertia forces, the deflections 
are obtained by an iterative procedure. For this purpose, 
a kineto-elastodynamic stretch rotation operator (EEDSEO) 
describing the elastic motions of a sii]g;le point is used. 

The crank is treated as a cantilever beam thus reducing the / 
mechanism to a structure and a system flexibility matrix is 
developed. To avoid excessive computations, Taylor aeries 
expansions of the matrices are used. These expansions 
effectively eliminate the deflections due to the elastic 
inertia forces. Erdman, Imam and Sandor'' ^ extended the 
above analysis by using the method of 'dynamic equivalence 
system' to include the elastic inertia forces . The authors 
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show that the elastic accelerations, in some cases, can be 
of the same order of magnitude as the rigid body accelerations, 

(17) 

Imam, Sandor and Kramer used the stiffness 
approach of the finite element method to analyse the planar 
flexible mechanisms. The assembly of the element matrices 
is done by permutation vector method generally reserved for 
large systems as it generates the system co-ordinates auto- 
matically within the program. It does not offer any special 
advantage because the system co-ordinates are less in number 
in case of mechanisms. On the other hand, it limits the 
flexibility of the program to some extent. The force vectors, 
calculated by multiplying the mass matrix with the negative 
of the acceleration vector, is assumed to remain constant 
during each interval. The errors in the force vectors, formed 
in this way, vary from 5% to 14,3% (in case of a purely 
rotating link) even at the element level, furthermore, these 
force vectors do not take into account the actual distribution 

of the rigid body inertia forces. The technique of the rate 

( 18 ) 

of change of eigenvalues and eigenvectors'’ ^ is employed 
for the first time in mechanism analysis and the computation 
time is reduced by a factor of three. However, it appears 

■m 

that the number of the eigenvalues to be taken in a mechanism 
without any intermediate check and the range of the crank ^ 
angle up to which the technique can be safely applied is yet 



10 


inconclusive in the general case. Mo system co-ordinate is 
defined corresponding to the element co-ordinate for rotation 
at the support-end of the crank. This effectively reduces 
the crank to a cantilever beam and consequently the mechanism 
to a structure. Though this procedure by-passes the 
necessity of eliminating the rigid body degree of freedom, 
it does not analyse the mechanism with its actual boundary 
conditions . 

1.1.2 Input torque analysis 

The fluctuation of the input torque due to the 

variable inertia forces have interested many researchers. 

But their investigations are ma,inly confined to the problem 

of balancing the ■ shaking forces and moments of the rigid 

(19 20 ) 

link mechanisms^ ’ . The importance of the fluctuation 

of the speed of a driving induction motor due to the variable 
moments exerted by a mechanism is first recognised by 
Habiger^^"' ^ , Muller^^^^ and Houben^^^”^^ ^ in their investi- 
gations of the torsional vibrations of the input shaft. 

It is shown that in some cases the fluctuation of the motor 
torque may enhance the instability of the system. 

In the areas of the planar rigid link mechanisms, 

( 2T PR 

Benedict and Tesar^ ^ found out the angular velocity 
fluctuations by employing the kinematic influence 
co-efficients and the following three methods; (i) energy 
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distribution method, (ii) equivalent mass and force method 

and (iii) rate of change of energy method. The fluctuating 

input torque can be next foimd out from the torque-speed 

characteristics of the driving motor, Erdman, Sandor and 

5 ) 

Oakberg extended this approach to the flexible mechanisms 
by including the rate of change of strain energy of the 
mechanism in the analysis. However, the accurate determin- 
ation of the rate of change of strain energy involves 

extensive computations. 

( 1112 ) 

Winfrey^ ’ has not incorporated any input 
torque in the force vector term of the equation of motion. 
Thus the unbalanced rigid body inertia forces must be 
balanced by the additional inertia forces due to the subse- 
quent elastic and rigid body accelerations. As the modal 
matrix does not include any rigid body mode and a constant 
input speed is assiuned, the subsequent rigid body acceler- 
ations during each interval are suppressed and the analysis 
does not represent the actual situation. 

Sadler and Sand or * assume that the input 
shaft has a large moment of inertia so that the fluctuation 
of the moment exerted by the mechanism on the input shaft 
does not affect the constancy of the crank speed and 

therefore the crank rotates as a cantilever beam. Imam and 
( 17 ) 

his colleagues also consider the crank as a cantilever 



beam as previously mentioned. This restriction on the crank 
obviates the necessity of finding out the variable input 
torque and its contribution towards the displacements of the 
mechanism. However, the moment at the cantilever-end may be 
easily found out by further calculations. 

1 .t .3 Harmonic analysis 

The importance of the harmonic analysis of 

mechanisms is borne out by the fact that it is the most 

efficient approach to describe the dynamic behaviour of 

mechanisms and to find out the characteristics of a motion 

or the shape of the path in a wide range of the mechanism 

parameters. It also serves as a powerful tool for the 

synthesis of mechanisms and enables the designer whether the 

designed mechanism possesses sizable overtones. However, 

the complexity involved in the harmonic analysis of the 

general mechanisms has been acting as an impediment to 

extensive research works in this area as evident from the 

scarcity of research papers published on this topic. The 

development of the harmonic analysis started with the planar 

slider-crank mechanism^ ^ and with the mechanisms involving 
(30-32) (33) 

slider pairs . Hreudenstein performed the harmonic 

analysis of a rigid link four-bar mechanism and proved that 
the satisfaction of certain requirements by the link lengths 
leads to optimisation of the transmission angle and 
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minimisation of the higher harmonics thereby reducing the 
accelerations, inertia forces, deflections and the dynamic 
bearing reactions. I'ollowing the publication of Freudenstein’ s 

work, a fev/ papers have come out on the harmonic analysis of 

/ *56) (37) 

planar mechanisms , spatial slider-crank mechanisms^ 

(38 39 ) 

and spherical four-bar mechanisms^ ’ > assiuming the 

links to be rigid. Sadler and Sandor^^^’'^"' ^ are the only 

authors who have attempted to extend the harmonic analysis 

in the area of the flexible mechanisms. The authors 

utilise the harmonic analysis of a crank-rocker 

(55) 

mechanism to find out the steady-state lateral vibrations 
of a concentrated mass attached to a flexible extension of 
the coupler. The crank and the rocker are assumed to be 
rigid; the coupler and its extension are considered massless 
elastic springs. Though the equation of motion for the 
longitudinal deflections of the mass is given, it is dropped 
out in the subsequent analysis by assuming the longitudinal 
inextensibility of the springs. The governing equation of 
motion is a non-homogene ous Hill equation in which the 
periodicity of the stiffness terms results from the changing 
geometry only. The unknown displacement is expressed as a 
Hourier series and the steady-state solution is found out 
by equating the co-efficients on both sides. Bending 
resonance (in both undamped and damped conditions) of the 
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concentrated mass is examined by numerical examples for a 

special case where the coupler is assixmed rigid. This work 

(41 ) 

is later extended^ to include the centrifugal force of 
the concentrated mass in the analysis. The analysis is very 
limited in the sense that the links are too much idealised 
and the stability analysis is not treated at all, 

1.1.4 Stress analysis 

When a continuum model is used to., study the 
deflections^^" f the normal and shear stresses within the 
links can be obtained without much difficulty from the 
stress-displacement relationships. Since all the authors 
using continuum model have concerned themselves with the 
deflections only, it remains to be seen how the stresses, 
involving the space -derivatives of the displacements, are 
affected when the actual displacement function is approxi- 
mated by one or two mode functions. 

(9 10 ) 

In the lumped-mass models of the mechanisms^ ’ , 

the stresses are calculated from the bending moments, 

axial forces and transverse forces, obtained in the course 

of deflection analysis. Since these procedures make use of 

the derivatives of the elastic curves, the stresses, thus 

(43) 

calculated, may not yield accurate results 
operating speed is high. 


when the 
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In the area of the finite element method applied 

to mechanisms, Beggs^^^^ first found out the stresses in a 

( 12 ) 

redundant mechanism for the static forces only. Winfrey 

does not find the stresses but only mentions about the 

computational difficulties in calculating the stresses from 

(17) 

the deflection analysis. Imam, Sandor and Kramer 
outline a new procedure for calculating the stresses within 
the links. This procedure assumes that when the mechanism 
vibrates in any of its normal modes, each link also vibrates 
(at least approximately) in one of its natural modes. This 
is not true in general. Moreover, the procedure requires 
that the deflections of the mechanism should be calculated 
and transformed back to the element co-ordinates separately ^ 
for each normal mode thus entailing complexities in the 
nume rical proce dure . 

1.1.5 Stability analysis 

Seevers and Yang^"^^’^^^ studied the stability of 

the elastic couple'-r of a slider-crank mechanism. The 

authors calculate the pin forces from the rigid body analys 

only and use only the first mode of deflection in the 

resulting Hill-type equations. 

( 47 — 49 ) 

Houben^ ' first investigated a four -bar 

mechanism treating all the moving links as elastic. The 
author develops nine simultaneous partial differential 
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equations with, periodic co-efficients which describe the 
bending in the mechanism plane as well as perpendicular to 
it and torsional vibrations but exclude the longitudinal 
deformations of the links. In the analysis, the author 
includes some linear terms coming from the relative tangential 
and normal components of accelerations, the Coriolis component 
of accelerations and the axial forces obtained from the 
rigid body analysis. Since the solution of these- simultaneous 
equations is extremely unwieldy, he assumes the driver and 
driven links to be rigid and, instead of solving the equations 
for the displacements, he studies the attendant instabilities 
by numerical methods. With the aid of an analog computer, 
the first stability zone is determined. 

Witfeld^^^^ investigated the stability and bending 
resonances of an elastic coupler of a planar four-bar 
mechanism for various combinations of linkage parameters. 
Calculating the pin forces from the rigid body analysis, the 
author finally forms a linear integro-diff erential equation 
with the curvature as the independent variable. Using a 
digital computer, he furnishes the stability and the resonance 
curves and discusses about the possibility of combination 
resonances . 

ful "5 

Tobias^ ■ studied the stability of an elastic 
coupler of a four-bar mechanism by giving small vibratory 
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motions to the crank about a rigid body position of the 
mechanism. The mechanism is idealised by considering the 
crank and the follower to be rigid. Moreover, the follower 
is assiimed to have a large moment of inertia compared to 
that of the coupler. Starting with Eringen’s nonlinear bending 

(55) 

equations^ , the author finds a power series expansion 

using the ratio of the inertia of the coupler to the inertia 

of the driven link as a parameter. The zero-th approximation 

is a four-bar mechanism with a massless but flexible coupler 

which becomes a spring-mass system with a time varying 

f 52 ) 

spring constant. Using Bolotin's results^ for parametric 
instabilities of a rod, the critical frequency zones for the 
axial instabilities are found out from the zero-th approxima- 
tion. Erom the zero-th and first order approximations, it is 
shown that the buckling can occur only when the axial load 
exceeds Euler's critical load. However, by perturbing the 
base solution, the critical frequency zones for the bending 
instabilities are found out. The analysis shows that (i) the 
axial vibrations of the coupler is of primairy importance and 
(ii) the bending instability occurs either when the axial 
load becomes sufficiently large to cause the column type 
buckling or due to the parametric resonance. 

It is proved that the bending instability is 
greatly influenced by the axial vibrations and the author 
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concludes that the axial vibrations should not be neglected. 

The similar conclusion has also been drawn in the studies of 

(■ 53 55 ) 

the parametric vibrations of columns and strings^ ’ . 

( 57 ) 

Broniarek and Sandor examined the stability of 
a parallelogram linkage with massless elastic coupler and 
rigid cranks, connecting two rotating masses mounted on elastic 
shafts. Because of these special properties of the links, 
the stability of the system is easily determined by studying 
the stability of a simple Hill-type equation. 

Some of the authors cited in Section 1.1.1 also 
discussed briefly in their papers the stability of the 
mechanisms. Mayer Zur Capellen^® ’ provides some resonance 
criteria and concludes that several natural modes can occur 
simultaneously for some input crank-speeds. Jainski, lee and 
Sandor^^’"^^ use the Routh -Hurwitz criterion to find out the 
stability zones of the coupler. Bor various frequency ratios, 
Viscomi and Ayre have plotted the transverse displacements 
where unbounded regions are shown but have not analysed 
further to study the stability problem. 

1.2 Objective of the Present Work 

Prom a careful study of the past research works 
cited in the previous section it is found that the solution 
with the continuum model is almost impossible except for 
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extremely idealised cases. Such models help only to provide 
an insight into the general nature of the problem. Several 
research works have been conducted to obtain a more comprehen- 
sive solution for more realistic cases by using numerical 
methods. The lumped mass model suffers from several drawbacks 
viz . , lack of generality in approach, complexity in the 
analysis procedure if the axial extensibilities are included 
and the inaccuracy in the calculated stresses. The above 
difficulties may be successfully eliminated by using the 
finite element model. However, several simplifying assumptions 
have been made in the previous works using the finite element 
models and till now no effort has been made to study the 
effects of the additional acceleration terms as well as the 
axial forces using the finite element models. It is evident 
that, of all the three models, the finite element model is the 
most powerful one and its potential can be exploited to obtain 
a comprehensive analysis of the mechanisms. 

The major objectives of the present work may be 
classified into three groups as follows s 

(i) The basic concepts of the finite element techniques 

will be utilised to provide a firm base for the analysis 
of mechanisms. The vibrating motions of a completely 
flexible mechanism will be studied by taking into 
account the axial vibrations of the links. To achieve 
this end, an effort will be made to examine the 



effects of subdividing the links on the accuracy of 
the solution. Unlike the previous works, the actual 
distribution of the rigid body inertia forces will be 
considered to generate the load vectors. It is also 
proposed to use a ‘ code -system’ of assembly to 
generalise the analysis which will be helpful in case 
of more complex situations. In addition, a general 
procedure for eliminating the rigid body degrees of 
freedom will be developed. As a natural extension of 
the above procedure, a simple equation will be 
evolved which gives the variable input forces of a 
general mechanism under certain assumption. With a 
modified assumption, this equation will be extended to 
study the vibrations of the flywheel at the crank- 
shaft. In this connection, a more practical case 
will be examined without imposing any simplifying 
assumption on the input forces. To study the dynamic 
behaviour of the rigid body inertia forces, the 
hamonic analysis of these forces will be made and 
their effects on the vibration characteristics of a 
mechanism will be investigated. A new method of stres 
analysis will be proposed which requires very little 
additional computations. It will also be demonstrated 
how the present procedure for the kineto-elastodynamic 
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analysis of a mechanism can be very efficiently 
utilised to obtain the rigid body analysis (viz., 
the kinematic analysis and the dynamic force analysis) 
of the mechanism. 

(ii) To assess the contributions of the various additional 
acceleration terms and the rigid body axial forces on 
the dynamics of the mechanisms, a general equation 
will be derived from Lagrange * s equation of motion by 
using a moving co-ordinate system. In order to avoid 
the complexity in the solution of the resulting 
equation of motion, an approximate method for solving 
these equations Virill also be outlined. An attempt 
will be made to study the effects of the elastic axial . 
forces on the displacements of the flexible mechanisms. 

(iii) In the conventional finite element method, the solutions 
are obtained only for an instantaneous configuration 

of the mechanism. Consequently, a large number of 
computations is needed to have a complete picture of 
the system behaviour. To overcome this difficulty, a 
new approach will be presented for a restricted class 
of mechanisms which will incorporate all the dynamic 
factors and describe the steady state displacements 
for all configurations of the mechanism by means of a 
single expression. As a result, the computation time 



for finding the steady state solution v\rill be reduced 
considerably. Since this form of solution is conducive 
to the stability analysis of mechanisms, the relevant 
procedures for determining the critical zones of the 
crank speed will also be described. 

The effectiveness of the above studies will be 
demonstrated by solving several problems and comparing the 
results with those from the past works whenever possible. 
However, due to the limitations of time available , the list 
of the numerically solved examples may not be exhaustive. 

For the same reason, no numerical example will be furnished 
corresponding to the following two procedures; for analysing 
the mechanisms with driving motors and for determining the 
instability zones of the mechanisms. 



CHAPTER II 


DEFIECTION AHD STRESS AHAiYSIS 

2.1 Introduction. 

It is apparent from the review of the past research 
that the published works on the finite element analysis of 
the flexible mechanisms lack the generality and rigour. In 
this chapter an attempt is made to place the basic analysis 
of the flexible mechanisms in a moderately rigorous and 
general framework of the finite element method. Section 2.2 
contains a precise statement of the problem and the assumptions 
made in the analysis. Section 2.5 outlines the analysis for 
each link at the element level. Section 2.4 describes a 
general method of obtaining the system matrices suitable for 
mechanisms. In Section 2.5 the solution of the system 
equations is carried out with a systematic approach to the 
problem followed by the determination of the deflections 
within the links. In Section 2.6 the problem of finding 
the input torque under various conditions is solved. A 
complete harmonic analysis of the forces for a slider-crank 
and a four-bar mechanism is performed in Section 2.7. A 
new and efficient procedure for the calculation of various 
types of stresses is presented in Section 2.8. The chapter 
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is concluded in Section 2.9 with a description of a new type 
of kinematic and force analysis of a rigid link mechanism. 

2.2 Statement of the Problem 

The problem to be dealt with specifically in the 
present study is to find out the deflections and stresses of 
a completely elastic mechanism under pre-ass igned external 
loading conditions. The problem also includes the exa- 
mination of a flexible mechanism for the limits of critical 
operating speeds. The major assumptions are: (i) the rigid 
body motion and the elastic vibrations are separable, (ii) 
the mean longitudinal vibrations of a link is independent 
of the transverse vibrations, (iii) the deflections are small 
so that linear theory can be applied throughout and (iv) 
the effects of the tolerances, clearances and impact are 
neglected. 

The equations and the method of analysis described 
in the present study are generally applicable for any 
mechanism. However, to make the presentation simple, a 
planar mechanism consisting of straight links with uniform 
cross-sections will be considered. 

To fit in the methodology of the finite element 
techniques, the following mathematical model of a mechanism 
is used. At an initial time x = x^, the rigid body 



configuration of the mechanism is considered as an 
’instantaneous structure’ capable of undergoing both the 
rigid body and elastic motions. (In case of a four-bar 
linkage, this position is shown by 0-|A^B^02 in Fig. 2.1). 

The forces acting on this instantaneous structure are the 
rigid body inertia forces due to the ’gross' rigid-body 
accelerations of the elements plus the external forces 
acting on it. The mass and stiffness properties of the 
mechanism, treated as an elastic structural system, are 
derived at this position and are assumed to remain unchanged 
during a chosen interval of time at. The displacements, 
stresses etc. for each element at the end of the interval 
are determined. From the deflected positions of the ends of 
the links, the shifted rigid body configuration (position 
is determined and is allowed to travel for a 
period of at as a rigid-link mechanism. The configuration 
thus reached (position 0 ^A^B^ 02 ) is considered as the 
instantaneous structure for the next interval of time at. 

The displacements and velocities of the elements at the end 
of this interval are taken as the ’initial conditions’ for 
the next interval. In this way, the entire process is 
repeated until the steady-state condition is reached. 

It is known that in dynamic problems, discreti- 
sation of a continuum gives accumte results only when 
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sufficient numbersof the element co-ordinates are defined. 

For this reason each link is assumed to be divided into 
several parts and each part is treated as an element of the 
mechanism. 

2.3 Element Analysis 

To start with, the element co-ordinate axes x-y 
are chosen along the axis of each element and perpendicvilar 
to it. Six element co-oix3inates U|('r) to ( t) are defined 
at the two ends of each element (denoted by E' ) along the 
corresponding element co-ordinate axes (Eig. 2.2) to represent 
the continuous behaviour of the element in a discretised 
way. The continuous axial displacement u^(x, y, -c) due to 
the axial forces only and the continuous transverse displace- 
ment Uy(x, y, t) due to the bending only can be approximated 
by the polynomials as expressed in eqs. (2.1 and (2.2). 

u^ = + CgX (2.1) 

Uy = *^3 + ^6^^ (2.2) 

where functions of time x only. 

The above polynomials are used because they 

^ Id 

relate the continuous displacements Uy to the nodal 
displacements u^ to Ug (Fig. 2.2) with the accuracy of the 
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engineering theory of static bending of beams. In dynamic 
problems other functions could as well be used and the choice 
of appropriate functions is a subject of research till now. 

It may be observed that eqs. (2.1) and (2.2) permit the 
rigid body motions in addition to the elastic motions. 

To include the effect of shear, the following 
relationships are used. 


9X GA 

O 

3 b 

El ^ = Stj, 


(2.3) 

(2.4) 


(2.5) 


where , 



the transverse deflection due to shear only 

total transverse displacement 

the shear force acting at a position x of 


the link 


A = the effective shear area of the cross-section 
s 

G = the shear modulus of elasticity 

E = the Young’s modulus of elasticity 

I = second moment of area of the cross-section 


Eliminating Sp from eqs, (2.4) and (2.5 ), the following 
equation is obtained? 
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( 2 . 6 ) 


Integrating the above equation and using eqs. (2.2) and (2.3), 
the following relation is obtained 


^ = o; . G,x . C^Cx^ + §1) . Cg(x5 . fl X) 

s s 


(2.7) 


Following boundary conditions are used to determine to 


Cg in terms of 112 to Ug 
At X = 0, u^^ = U2J 


7 


At 


X = 1, 


= u. 


Uy - 


8 x - ^3 


ii-5 

ax ~ ^6 


(2.8) 


Finally, 

«y = (1 "I ”) ■ [l-3g^+2g^+(l>g)9 l-2?2+5^4i(?-5^)9 


3 e^- 2 C^+ 9 C -C^+5^-f(?-5^)] 


Ur 


u^l 


U^ 




(2.9) 


1 2EI X 

where, 9= 5 5 = r- > 1 = length of the element. 

GA^i2 

Id 

The longitudinal displacements due to the 
bending is given, in accordance with the engineering bending 
theory, by eq. (2,10), 
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uj = - ^ y = {-i+ 4 { -352.(1. c^jc 


-(65-65^)s (26 -35^-?«)c] 


u. 

c. 

Ujl (2.10) 

^5 

L^e^J 

With 5 = y /1 

Using the boundary conditions u^(0,t) = u^, -0^(1, t) = u^, 

and are determined from eq, (2.1). The total displacement 
Ujj. in the x direction is 

u^(x,y,T) = u|(x,t) + u^(x,y,T) (2.11) 

Thus, the continuous displacement ^ ^y) expressed 

in terms of the nodal displacements u = {u^ , ..., Ug} in eq. (2.12a) 

u=au (2.12a) 

Consequently, u=au; u=au (2.12b) 

where the shape function a is given by eq. (2.13). 

1 

TV I 

0 1-35^+2C^+(1-c) 


a ?s: 




(1+(j>)(1-5) 6(?-c‘^) X, 


U,-l+ 4 g- 552 _('i.. ^ (l+(f>)5 6(-?+r) ? (25-35 -Cf)lC 


[ 5 -25^+5^+i(5-5^)« JL 


0 


55*^-25^+ [-5^+5^“2 (c~€ ) ^ P- J_ 


_2x 


(2.13) 




with c = y/i* 


The normal strain e^„Cx, y, x ) and the shearing 

strain e (x, y, x) at any point of the link E' are given by 
xy 

;^he following equations: 
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(2.U) 

(2.15) 


Introduction of eq. (2.12) in eqs. (2.14), (2.15) gives the 
continuous element strains e ={e , e Iwithin the element 

E' in terms of u as shown below: 

e=bu (2.16) 


where the element strain matrix b is 
b. 


1 


(ITtJI 


-(1 + f) 6(1~2c)? (4-65+9)15 (i+f) 6 (-1+25) 5 (2-65-9)1^ 
0 -9 - ^ 0 9 


“! 


_ 

2 


(2.17) 


Expressions for normal stress h^ 5 j-(x, y, v) and nominal shear 


stress a (x,y, t) at any point of the element E' are given 
xy 

by eqs. (2.18) and (2.19) respectively. 


a = E e 

XX XX 


(2.18) 


xy 
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I 


(2.19) 
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Using eq. (2.12) in eqs. (2.16) to (2.19). continuous 
element stresses a = within the element E’ are 

obtained in terms of u as follows? 


where 



( 2 . 20 ) 


( 2 . 21 ) 


The kinematic analysis for a rigid link mechanism 
can be performed either by iteratiwe procedures ^ ^ ^ or by 

the method described in Section 2.9. For illustrative 


purposes, direct kinematic analysis for two planar mechanisms 
is given in Appendix D, 


From the kinematic analysis, the absolute rigid 
body normal accelerations A^^, and tangential accelerations 

Aj^y, A^^ at the end points M and U respectively of the element 
E' are determined. The components of the corresponding body 
forces B(5, x), varying linearly within the ‘element, are 
given by eq. (2.22). 



where p = density of the material. 



M 


( 58 ) 

I’rom the principle of virtual work, the element 
oriented element stiffness matrix k, the element mass-matrix 
m and the element load -vector ^ are obtained as shown below; 


"E = / b"*" X b dv 

V 

ffi = / pa"^ a dv 

V 

p = / a"^ B dv 

V 


( 2 . 23 ) 


( 2 . 24 ) 


- / p a^ 

V 


"Sax ^ 


( 2 . 25 ) 


where t denotes the transpose of a matrix and v implies 
that the integrations are carried over the whole volume of 
the element . 

Derivation of these matrices and their explicit 
expressions are included in Appendix B. 

To define the properties of the element E' in a 
common co-ordinate system, a datum (or system or global) 
co-ordinate system X-Y is chosen and a rotation matrix R is 
used for the transformation of the element oriented element 
co-ordinates into the system oriented el men t co-ordinates. 
The matrix R-'can be expressed as follows: 

" R.J ; 0 - 

0~t^ 


R 


( 2 . 26 ) 
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where is the rotation matrix for each end and is given 

by the following equation (for plane mechanisms): 

“ cos©„, sin0_, 0” 

E' E* 

R^ = - sinG-^, cos©^, 0 (2.27) 

0 0 1 _ 

By virtue of the oo-ordinate transformation defined by eq. 

-IS. A 

(2.28), the load vector p, the stiffness matrix K and the 
mass matrix in the system oriented element co-ordinates 
are obtained, from the contragradient law of transformation, 
in eqs. (2.29) to (2.31 )• 


u = 

R u 

(2.28) 

.'A 

k = 

r"^ k R 

(2.29) 

m = 

R*^ m R 

(2.30) 

p = 

r'^ p 

(2.31) 


where u is the nodal displacement, vector in the system 
oriented element co-ordinates. The initial displacement 
u(t.^) = u*^ and initial velocity S(t^) = u in the system 
oriented element co-ordinates are obtained in eqs. (2.32) 
and (2.33 ) from the corresponding terms u(t:q) =13- and u(Tq) = u° 
in the element oriented element co-ordinates. 




^0 

u 

= r'^ 

if 

(2.32) 


= R^ 


(2.33) 


-fj 

since R R = I, the identity matrix. 

If there is any externa.1 distributed or concentrated 

force acting on the element, it must be converted into an 

equivalent element force vector p^, derived from the equivalence 

of virtual work^^^\ Ror example, if a surface load p acts 

s 

on the element over a length extending from x = s^ to x = s^, 
the equivalent element force vector p is given by eq. (2.33a). 

Pg = J a'^ Pg ds • (2.33a) 


In case of the concentrated forces p^, s^ in 
eq. (2.33''a) approaches to s^ and p„ds becomes equal to p_ in 

o V 

the limit . 

The equivalent force vector p found out- by 

,j frj 

eq. (2,33a) is conveniently added either with p or with p 
depending on its direction. 

The above element analysis is carried out for all 
the elements present in the mechanism and the element 
analysis is completed. 
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2.4 Assembly of Element Matrices 

The assembly of the element matrices (found out 
in the last section) into system matrices for the whole 
mechanism is carried out by 'code-system' method^ . 

To explain the method, identification numbers 
(encircled in Eig. 2.3) are alloted to each element of the 
mechanism at the outset. An integer matrix US is defined 
such that its columns contain the system co-ordinates 
associated with the system-oriented element co-ordinates of 
an element whose identification number matches with the row 
number of these columns. For example, US for the mechanism 
shown in Fig, 2,3a is given by eq. (2.34). 

11 12 10 1 2 9 
US:= 1 2 3 4 5 6 

5 6 7 13 14 8 

let the element number 1 in Fig, 2.3a be assumed 
rigid but allowed to rotate while in Fig. 2.3b, the element 
number 1 be assumed not only rigid but temporarily fixed in 
space during each interval. For that ptxrpose, the total 
number of the degrees of freedom for these mechanisms are . 
specified as 8 and 6 respectively and the. system co-ordinates 
are numbered accordingly as shown. When no information like 
the reactions etc. at the co-ordinates with zero displace- 
ments (Fig. 2.3c and Fig. 2.3d, for illustrative purposes) are 


(2.34) 








required, these co-ordinates can he alloted the same number 
to reduce the size of the matrices. 

In this method of assembly, the number of the 
element co-ordinates and the system co-ordinates at any 
type of joints need not follow a regular order. For example, 
joint G (Fig. 2.3c) having three branches and element 2 
(Fig. 2.3a) having 3 nodes but no moment co-ordinate- can 
be handled without any difficulty. Moreover, the dependent 
co-ordinates (defined afterwards) can be grouped together 
as the last set among the generalised system co-ordinates 
and the advantage of this type of flexibility in numbering 
the co-ordinates is manifested in a later section. As the 
association of the system co-ordinates with any element can 
be readily established by simply referring to NS , the special 
features of a particular mechanism, viz., the mass and the 
centrifugal effect of the concentrated weight at F in Fig. 
2.3<3 .j the gas force, the piston mass and its inertia in 
Fig. 2.3b, the input and output torques etc. can be directly 
incorporated to the system even when the division of the 
links are made arbitrarily. The basic difference between 
the conventional method^^^’ and this method of assembly 
is that in the former case the system co-ordinates are 
numbered automatically inside the program thus imposing 
certain restrictions in its use while in the later case it 
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is done through the input data rendering complete freedom 
in its application. Unlike other structures, the total 
number of the elements and joints in a mechanism is small; 
so the automatic numbering does not offer any special 
advantage over the little disadvantage of using a relatively 
bigger matrix US. On the contrary, it is rather advantageous 
to feed the boundary conditions, regarding the types of the 
joints, in the program directly through a single matrix IB . 

In view of the above, the code-system seems to be more 
general and efficient compared to other methods of assembly 
like connectivity matrix \ permutation vector method 
etc. used in the earlier works. 

A sample program to perform the assembly in case 
of a flexible mechanism is given in Table 2.1 where the 
links can be subdivided arbitrarily. 

Using the code-system as described above , the 
system, matrices are generated from the corresponding matrices 
in the system oriented element co-ordinates as indicated 
in Table 2.1. Having done this, the special features of the 
mechanism under consideration are incorporated in the system 
matrices, for example, the inertia of the flywheel at the 
rocker end (fig. 2.3'd) should be added to a diagonal element 
of the system mass matrix. The row (or column) number of 
this diagonal element corresponds to the number of the 
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Table 2.1: Sample program for assembly 


IDS = 0 
DO 1 I = 1,M 

IDIV = DIV(I) 

CALL ESEM(I,LS,EM,SP) 

DO 1 miY = IDIV 

BIDS = ms + 1 

DO 1 J = IjMEC 
lEOYf = IS(MDS,J) 

SP(IROW) = SP(IROW) + EP(J) 

SUO(IROW) = EUO(J) 

SVO(IROW) = EVO(J) 

DO 1 K = 

ICOL = RS(lviDS,K) 

SS(IROW,ICOL) = SS(IROW,ICOL) + ES(J,K) 

1 SM(IROW,ICOL) = SM(IR0W, ICOL) + EM(J,K) 


M = Riimber of linksi DIV(I) = number of divisions for the 
i-th link; mC = maximum number of element co-ordinates used. 
ES, EM, EP, EUO, EVO and SS, SM, SP, SUO, SVO are ssrstem 
oriented element matrices and system matrices for stiffness, 
mass, load, initial displacement and initial velocity 
respectively. It is assumed that the elements of a link are 
numbered sequentially. 
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system co-ordinate associated with the rotation of that 
rocker end. Siinila.rly, the centrifugal force and the dead 
weight of the concentrated mass at F (Fig. 2.3d) should he 
added to the elements of the system's force vector whose 
row n-umhers correspond to the numbers of the horizontal and 
vertical system co-ordinates respectively at the point F, 
while its mass should be added to the two corresponding 
diagonal elements of the system mass matrix. Also the 
forces or torques exerted by the input shafts (determination 
of these input forces is discussed in Section 2.5.3) are 
added to the corresponding elements of the force vector of 
the system. 

When the special features of the mechanism are 
incorporated, the final equation of motion for the whole 
mechanism is given by eqs. (2.35) and (2.35a) 'respectively. 

MU+CU+KU=P ■ (2.35) 

for -r^ < T < v Q + A'c ♦ 
subject to the initial conditions: 

U(t^) = and U(t^) = (2.35a) 

• 

where M, C, K, P, U(r), U° and U° are the mass matrix, 
damping matrix, stiffness matrix, force-vector, displacement 
vector, initial displacement vector and initial velocity 
vector respectively for the whole mechanism expressed in 



the system co-ordinates. Determination of the damping matrix: 

C is deferred to Section 2,5.4. 

2.5 System Analysis 

2,5.1 Derivation of the transformation matrix 

To solve eq. (2.35) by normal mode method for the 
first few modes, it is necessary to solve the corresponding 
eigenvalue problem as defined by eq. (2.36) 

M X + K X = 0 . (2.36) 

The solution of eq. (2.36) by Power method 
requires the inversion of K. Por mechanisms, the stiffness 
matrix K is singular^^^^ and the degree of singularity (or 
nullity) of the matrix: K is equal to the number of the 
rigid body degrees of freedom present in the mechanism. A 
general method to render K non-singular by removing its 
rigid body degrees of freedom is presented below. 

Prom the total set of the co-ordinates X of eq. 

(2. 36),, any subset of the dependent co-ordinates Xg is 
chosen in such a way that no rigid body displacement in the 
mechanism is possible through the remaining independent 
co-ordinates Xj^, It is clear that the total number of the 
co-ordinates present in X^ must be equal to the number of 
the rigid body degrees of freedom of the mechanism and 
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many combinations are possible to form the set of co-ordinates 
X^. However, it is convenient to choose the co-ordinates 
corresponding to the input forces as dependent co-ordinates 
and the advantage of it will be apparent in Section 2.5.3. 
Accordingly eg.. (2.36) is rewritten in partitioned form as 
shown below. 











%A 



« • 


%A 

^AB 


^A 

J 

..J 


1 


Kba 

%b 


^B 


0 (2.37) 


If any arbitrary small rigid body displacement 
dXg is given to the dependent co-ordinates, the resulting 
rigid body displacements dX^ at the independent co-ordinates 
may be expressed in terms of dX^ in the following form, 

dX_^ = I dX^ (2.38) 


I in eq,. (2.38) is called the displacement 
transformation matrix*"^^^ and depends only on the geometry 
of the mechanism. Ihe procedure of determining I directly 
from the geometry in case of simple mechanisms is exemplified 

in Appendix 0. A general procedure for the determination 

■ 

of T in case of any mechanism is outlined next. 

Since a mechanism does not produce any elastic 
force due to the rigid body displacements, e<ls. (2.39) and 
(2.40) immediately follow. 
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That is , 


^AA ^XA ^AB ^XB “ ^ 

(2.40) 

^A ^XA ^B ^^B = ° 

Introduction of eq. (2.38) into eq. (2.40) yields 
the following useful relations; 



Equa.tion (2.41) offers a most general method for 
finding T and eq. (2.42) can be used for cchecking purposes. 

It is interesting to note that though material and cross- 
sectional properties are involved in and K^-g, the final 
expression of T solely depends on the geometry of the 
mechanism. Eq. (2.41) is most efficiently exploited if the 
inversion of is utilised in the solution of the eigenvalue 
problem defined by eq. (2.36). Eor this purpose, a new 
procedure is described below which requires the inversion 
of only. 



2.5.2 Elimination of the rigid body degrees of freedom 


Let any virtual rigid body displacements be 
given at the dependent co-ordinates of the system defined 
by the equilibrium equation (2.37). fhe resulting virtual 
rigid body displacements 6X^ at the independent co-ordinates 
follow from eq. (2.38) and is expressed in eq. (2.43). 


6Xa =1 6X, 


(2.45) 


Since the restoring forces do not produce any work 
on the rigid body displacements, work is done by the inertia 
forces only. Thus, the following equation is obtained from 
the principle of virtual work when applied to the equilibrium 
equation (2.37). 

+ %B*b) + + Vb' = ° <2.44) 

Use of eq. (2.43) and the fact that SX^ is arbitrary, yield 
eq. (2.45) relating Xg with X^. 

Xg = T'* X_^ (2.45) 

where T* = - + Mgg) ^ (2.46) 

/ CO Y 

For small displacements, can be resolved into two parts^ ^ 
as shown in eq. (2.47). 


( 2 . 46 ) 




( 2 . 47 ) 


47 


In eq. (2.47) denotes the displacements at the 
independent co-ordinates relative to the rigid body configu- 
ration corresponding to Xg = 0. Direct use of eq. (2.45) in 
eq. (2.47) to eliminate X^ yields the following equation 
relating X_^ with X^^. 

X^ = (I - T T*)"^ X^^ (2.48) 

The matrix I - T T* is usually large and therefore its 
inversion is time consuming. Alternatively, eq. (2.45) can 
be used in eq. (2.47) to relate X^ with X^^, after the 
elimination of X^, resulting in the following equation. 

= (I - T*T)~^ (2.49) 

The size of the matrix I - T*T is equal to the total number 
of the rigid body degrees of freedom and its inversion is 
much simpler than that of I - T T* . 

is eliminated from eq. (2.37) with the help of 
eq. (2.47) and from the resulting equation, X^ is eliminated 
by using eqs. (2.41) and (2.49) to form the following eigen- 
value problem expressed in the co-ordinates of X^^^. 

+ ^AA^^Ae = ^ (2.50) 


where , 



iS 


CO = natural frequency of vibration of the system 
M‘ = effective mass matrix in the free-vihration 
of the mechanism 

= “aA + ^“aa^ + ^ ' (2.51) 

If n and w be the total number of the system co-ordinates 
and the rigid body degrees of freedom of the mechanism 
respectively, the size of the matrices ^-ad M' is 
(n-w)X(n-w) . 

For single degree of freedom systems, I - is 

scalar and M’ is symmetric. Since is non-singular, eq. 

(2.50) is solved by Power method for the first few (say 
m in number) natural frequencies and natural modes. The 
modal matrix 0’ corresponding to eq. (2.50) is thus found out. 
Using eqs . (2.47) and (2.49), the modal matrix ^ correspon- 
ding to eq. (2.37) is found out as follows; 

[^-] (2-52) 

B 

where , 

= (I - T*T)~’’ (2.55) 

and = $’ + T (2.54) 

Besides the utilisation pf ihe above 

eigenvalue formulation in the co-ordinates of requires 
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lesser amount of computations compared to similar formulations 
with respect to or in a sub-space of and ^ . 

2.5.3 Variable input forces 

Since the inertia forces arising from the gross 
rigid body moTements of the mechanism vary at every position, 
the input shafts must exert variable reactive forces (which 
may include torques also) to maintain the stipulated input 
conditions at different positions of the mechanism. The 
problem of finding these variable input forces in a practical 
operating condition. is examined in Section 2.6.2. for the 
time being, it is assumed that the input shafts exert forces 
Tg at the dependent co-ordinates of the mechanism in such a 
way that the rigid body inertia forces together with the 
external forces acting on the ’instantaneous structure’ are 
in 'static' equilibrium in the rigid body degrees of freedom. 

If and fg are the total forces (the rigid body 
inertia forces plus other external forces except Tg) acting 
at the independent and dependent co-ordinates of the mechanism 
respectively, the corresponding force-displacement relationship 
can be expressed by eq. (2.55) in partitioned form. 


^AA 



1 


^A 

1 

Sb _ 


1 


^Bg + Tg 


(2.55) 



If any virtual rigid body displacements 6Xg be given at the 
dependent co-ordinates, the resulting virtual rigid body 
displacements at the independent co-ordinates are obtained 
as before from eq. (2.43). A,s mentioned earlier, the elastic 
forces do not produce any work in this case and eq, (2.56) 
follows from the principle of virtual work. 

6X^%A + ® (2.56) 

Eliminating dX^ from eqs, (2.45) and (2.56) and noting that 
6Xg is arbitrary, the following equation is obtained: 

Tg + Eg + T^Ea =0 (2.51) 

Eq. ( 2 . 57 ) expresses the equilibrium equations 
(total mjrnber of these equations being equal to the total 
number of the rigid body degrees of freedom) which are 
satisfied by the self-equilibriating forces acting on the 
mechanism. 

Input forces, thus found out from eq. (2.57), are 
added to the corresponding elements of P in eq. (2.35). 

2 . 5.4 Solution of equations 

Though the displacement vector U in eq. (2.35) 
is capable of executing the rigid-body motions as well as^^^^ ^ ^ 
the elastic motions (by virtue of the polynomials in eqs . 

(2.1) and (2.2)), the modal matrix $ contains only the 
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elastic modes and no rigid body modes. However, it was 
assumed in Section 2,5.3 that the rigid body forces including 
the input forces acting on the ’instantaneous structure* 
balance among themselves. Under this assimption, it is 
shown in Section 2,6*2 that no rigid body motion is possible 
for the ’instantaneous structure’ during the interval At. 
Therefore the moda.1 matrix $ , as determined in eq. (2.52), 
is used for the solution of eq. (2,35). A case where it is 
necessary to include the rigid body motions in $ is examined 
in Section 2,6.2. 

To solve the equation of motion (2.35) by the 
normal mode method, the system co-ordinates of U are trans- 
formed to m numbers of normal co-ordinates as shown by eq. 
(2.58a). 

U = 2>'n (2.58a) 

where p(t) is the displacement vector in the normal co- 
ordinate system. 

Noting that 5 is constant during the interval At, 
successive differentiations of eq, (2.58a) with respect to 
time T yield 



U = 

®f] 

(2.58b) 

and 

ij = 


(2.58g) 


l.i.T. r/'f'-FUR 
CEHTKAL LIBRARY 


?fhere f) and Yi are the velocity and acceleration vectors in 
the normal co-ordinate system. 

Substituting eqs. (2.58) in eq. (2.35) and pre- 
multiplying it by eq. (2.59) is obtained. 


Mt) + Of) + Ktl 


(2.59) 


where , 


M = 


K = ^ K $ 




P = 


( 2 . 60 ) 


Employing eqs. (2,58), the initial conditions for eq. (2.59) 
are expressed in the following form. 


(jh)-'' 

u° 

(2.61 ) 



(2.62) 


Because of the orthogonality of the modes in $ , M and K are 
diagonal matrices. Obviously, the matrix 0 is also diagonal 
if the damping matrix C is taken as proportional to either 
M or K or a linear combination of them. But the use of a 
single constant (proportionality constant) to represent the 
damping characteristics of the systems with multiple degrees 
of freedom is not realistic. Por this reason, the standard 
practice is to consider C as proportional to the critical 
damping of the system at the normal modes . ^ 


For this 
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case also, G becomes diagonal and equation- (2.59) become 
uncoupled. These uncoupled equations assume the following 
form 



where , 


= displacement at the i-th normal co-ordinate 
f|^ = velocity at the i-th normal co-ordinate 

= acceleration at the i-th normal co-ordinate 
m^ = i-th diagonal of M 
= i-th diagonal of K 


p . = i-th element of P 

= i-th natural frequency 


2m. ea. = critical damping at the i-th normal mode 

= damping ratio at the i-th 


normal mode 


As the damping present in a system tends to 
suppress the higher modes, it is desirable to take the higher 
values of Cj^for the successively higher modes. 

Por the choice of the damping shown in eq. (2.63), 
the damping matrix C is given by eq. (2.64). 




c 



Hi K 

m m 


M 


(2.64) 


Therefore, the system damping matrix G is eventually obtained 
from eq, (2.60) as 




(2.65) 


As the damping is chosen to be proportional to the critical 
damping, the full f oms of the damping matrices C and C as 
shown by eqs. (2.64) and (2.65) are never required in actual 
computations . 

Section 2.7 describes the solution procedure wten 
P ineq. (2.35) varies with time. However, in this Section 
P (and therefore P) is assumed to remain constant during the 
interval At • Por this case, the solution of eq. (2.63) is 
given by eq. (2.66). 


Ti^(T) = e 


— or. X ^ p . 

^ (I. cospcT’ + J- sinjj,. t*) + — ^=“2 

^ ^ ^ m.(D. 

XI 


( 2 . 66 ) 


Consequently, 


•cr. x ' 


f) (v) = e ^ + p^Jj^)cosp^TV - (o^J^ + p^Ij_)sin 


(2.67) 






(2.69) 


T)? = i-th element of the vector ti° in eq. (2.61) 
q ^ = i-th element of the vector r)° in eq. (2.62) 

The displacement p, velocity f) and acceleration 1), 
determined from eqs, (2.66) to (2.68), are transformed hack 
first to the system co-ordinates with the aid of eq. (2.58) 
and then to the system oriented element co-ordinates for each 
element by using the integer matrix US (as discussed in 
Section 2.4). These are finally transformed to the element 
oriented element co-ordinates by using the transformation 
given by eq. (2.28). The accelerations, transformed in this 
way, are used in Section 2.8 for the stress analysis 
of the mechanism. 



It may be noted that the element oriented element 
displacements + A't) and element velocities u(t^ + at), 

thus found out, are measured in the element oriented element 
co-ordinates and they are the only co-ordinate systems which 
maintain a time -invariant relation with the corresponding 
rigid elements. For this reason, u and u, as determined 
above, for each element are used as the 'initial conditions* 
u° and u in the eqs. (2.32) and (2.. 33) respectively for the 
next interval of time at, Next, starting from eq. (2.26), 
the whole procedure is repeated for the next interval after 
updating the initial time t^ 

2.6 Input Forces for Different Conditions 

2.6.1 Input forces for 'dynamic' equilibrium 

In this section, the analysis described in Section 

2.5.2 is modified by making the following assumption; the 
variable input forces balance the rigid body forces 
(inertia and externally applied) together with all elastic 
inertia forces acting on the instantaneous structure, i.e. , 
the instantaneous structure remains in 'dynamic' equilibrium 
for the rigid body degrees of freedom. 

The equation of motion (2.35) for the whole 
mechanism is rewritten in partitioned form in eq. (2,70) 
where T-o Is shown separately as an unknown quantity in the 
force vector P. 
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^BA Sb J L^B 




(2.70) 


Since the mechanism acted upon hy the inertia 
forces -mi, Tiscous forces -CU and rigid body forces P m 
eq. (2.70) is assumed to be in equilibrium, these for 

must satisfy eq. (2.57). Hence, 

- Va - Wb - '=BA°A - Wb + ^-“aa'^A 

- '''ab'’b - °AA°A - °AB°B = ° 


(2.71) 


or, l-t 


- - !b + ^“ba + ’^^aa^'^a + ^“bb + 


+ (°BA + ®^°AA^^A ^“bB ® “aB^^B 


(2.72) 


Upon substituting eq. (2.72) in eq. (2.70), the 
following equation of motion in modified form is obtained^ 




-t'^c 


-I'^C 


Kba ^ab 


%A ^BB J pB 


-t'^p 


(2.75' 
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It may be noted that in this case also the special 
form of the force vector in the right side of eq. (2.73) 
obviates the necessity (explained in Section 2,6.2) to 
include the rigid body inodes in the solution of eq. (2.73), 
Therefore eq. (2.73) can be solved in a manner similar to 
that used for e^. (2.35) in the Sections 2.5.3 and 2.5,4-, 

2.6.2 Input forces supplied by motors 

In this section a more practical case is considered 
where the input forces are given beforehand (and therefore 
not to be determined from the equilibrium considerations). 
Though the procedure presented here is quite general, only 
a particular case is treated for the clarity of presentation, 
let the power to drive the mechanism be supplied 
by induo t ion motors attached to the massless rigid input 
shafts. T’or the known speeds of the input shafts, the 
torques exerted by the motors at that instant are known from 
their torque-speed characteristics. As these torques in 
general do not balance the total moment exerted by the rigid 
body forces acting on the instantaneous structure, the 
rigid body motions (in addition to the elastic motions) take 
place during the interval Ax caused by the net unbalaneed 
moments. Thus, it is necessary in this case to include the 
rigid body modes in the modal matrix ® in eq. (2,.58) to alloxv 



59 


the ::rl.gid body motions during the solution of equation of 
moiii-om (2.35). 

It is assumed that the torques exerted by the 
mo toaxs (determined from the speeds of the input shafts at 
"tbie heeginning of the present interval At) remain constant 
duarilngg the interval and are added to the elements in the 
apprro;priate rows of the system force vector P in eq. (2,55). 
AILS'O it is assumed that the polar moments of inertia of the 
rot orrs are added to the appropriate diagonal elements 
( coarrresponding to the rotational co-ordinates of the input 
e ndls of the cranks) of the system mass matrix. 

It is clear that the modal matrix ¥ corresponding 
to tFlie rigid body modes _ of the instantaneous structure is 
gi'nreia by the following relation: 



(2.74) 


Replacing {30^, } and { 6X^, 6Xg} in eq. (2.44) 

bjr M and ¥ respectively, equation (2.44) is rewritten in the 
fonm shown below: 

M $ = 0 (2.75) 

Noting that the elastic forces do not produce any 
w orrk on the rigid body modes, eq. (2.76) follows. 


t y 

¥ K $ 


0 


( 2 . 76 ) 



Eqs . (2.75) and (2.76) express the orthogonality 
conditions between ^ and ^ with respect to the mass matrix M 
and stiffness matrix K respectively. Since dXj^ , dX^ in 

eq. ( 2 . 39 ) is any rigid body mode, the following equation is 

# 

obtained; 

K ¥ = 0 (2.77) 

In general, the rigid body modes in are not 
orthogonal to each other with respect to the mass matrix, 

i.e . , 

M f = M° (2.78) 

where is a w X w matrix (non-diagonal in general). 

In contrast to the relations among the elastic 
eigenmodes (i.e., they are orthogonal to each other with 
respect to both the stiffness and mass matrix), each of the 
rigid body modes are orthogonal to each row of the stiffness 
matrix, as seen from eq. (2.77), while they are not in 
general orthogonal to each other with respect to the mass 
matrix as shown by eq. (2.78). 

Recalling the special advantage of the normal 
mode method of analysis (i.e., the final solution being the 
superposition of the individual solutions obtained for each 
mode), the solution of eq. (2.35) is shown below for the 
rigid body modes only as the solution for elastic modes 
has already been described in Section 2.5.4. 
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Let w numbers of the rigid body displaoements in 
the normal co-ordinates corresponding to w numbers of 
the rigid body modes be denoted by the vector xix ) . The 
following transformation of the co-ordinates is used to 
describe the rigid body modes of motion by the normal 
co-ordinates. 

U = YX (2.79) 

Upon substitution of eq. (2.79) in eq. (2.35) and- premulti- 
'fc 

plying it by f , the following equation is obtained. 

M f X + K X = P (2.80) 

As the material damping does not occur in the rigid body 
modes, the damping term drops out from eq. (2.80). 

By virtue of eqs, (2.77) and (2.78), eq. (2.80) 

reduces to 

X = (M®)""' P (2.81) 

Assiming that P is independent of time, the solution of eq. 
(2.81) is easily obtained as follows: 

X(t) = ^(M°)”^ P + At + B (2.82) 

where the constants of integration A and B are to be deter- 
mined from the initial displacement vector x(Tq) and 
initial velocity vector x(Tq). 



The concept of the instantaneous structure is 
introduced in the present study (explained in Section 2.2) 
to perform a step by step analysis of the mechanism by 
dividing its continuous motion into numerous successive 
intervals of time. At the beginning of each interval, the 
total motion, inherited from the previous time histoiy, is 
separated into two parts: (i) the rigid body motions 
(represented in the form of the rigid body inertia forces) 
constituting a part of the force vector P in eq. (2.35) and 
(ii) the elastic motions represented by the initial elastic 
displacement U° and the initial velocity U° constituting the 
'initial conditions' of eq. (2.55a). The current rigid body 
configuration of the mechanism (it is rotating with respect 
to the fixed frame of reference) is taken as the datum 
position and during the interval under consideration., the 
mechanism is allowed to execute the rigid body motions as 
well as the elastic motions, both measured from this datum 
position. Thus the rigid body components of the initial 
displacement U° in eq. (2.35a) and consequently ^(t^q) is zero. 
Let the rigid body velocities in the system co-ordinates, 
determined from the kinematic analysis, be denoted by the 
vector Then using a transformation analogous to eq. 

(2.79), the initial velocity vector X (t^) = X is given by 



(2.84) 


After introducing these initial conditions, eq. 

(2.82) reduces to; 

X(t) = P ^ *,^0^ 

and X(t) = (M°)""' p r + X° (2.84a) 

By back substitution of eq. (2,84) into eq. (2.79), the rigid 
body displacements at the system co-ordinates at the end of the 
interval hx is found out. Similarly, the rigid body velocities 
at the system co-ordinates are obtained by back substitution 
of eq. (2.84a) in an equation obtained by differentiating 
eq. (2.79) with respect to time. 

These rigid body displacements, velocities given 
by eqs. (2.84) and (2.84a) are finally transformed into the 
element oriented element co-ordinates following the same 
procedure as described in Section 2.5.4, i.e., by integer 
matrix NS and eq. (2.28). The rigid body displacements and 
velocities so obtained at the input shafts are next added to 
the ’gross’ rigid body displacements and speeds of the input 
shafts prevailing at the beginning of the interval. The 
total velocities of the input shafts, thus determined, are 
used to find out the input torques for the next interval. 
Similarly the total rigid body displacements at the input 
shafts, as determined above, are used for the rigid body 
kinematic analysis for the next interval. 



In some positions of the me-ohanism, the 
rigid body moments exerted by the mechanism may differ 
considerably from the torques exerted by the motors and may, 
in turn, produce relatively large rigid body accelerations 
at those positions. It may be noted, however, that the 
theory of the finite elements described here is valid for 
small displacements only and, therefore, proper care should 
be taken to choose at sufficiently small at these positions 
so that the rigid body motions produced in these intervals 
remain within limits. 

As X in eq . (2.84) represents only the small rigid 
body motions with respect to rigid frame of reference, the 
’gross’ rigid body displacements must be found out each time 
separately from the rigid body kinematic analysis. The 
rigid body velocities can be merged with elastic velocities 
U° to form a total initial velocity vector. Because of the 
orthogonality relations between the rigid body modes and the 
elastic modes with respect to the mass matrix M, the 
respective components of this total initial velocity vector 
can be isolated by the matrix operations''^ . However, the 
modes, determined by the numerical procedures, are orthogonal 
to each other only within a certain degree of accuracy. 
Therefore the isolation process will involve some inaccuracy. 
This inaccuracy may accumulate in the successive intervals 



and may cause substantial error in the final results. 

Moreover, the rigid body velocities may vary considerably from 
the elastic velocities and even slight inaccuracy in the 
modes will affect the separation process adversely. Therefore 
it is better to treat separately from u'^ as described 
above, ■ 

It is seen that if the force vector P is self- 
equilibriating, the first terms in the right side of eqs. (2.81), 
(2.84) and (2,84a) become zero by virtue of eq. (2.74). 

Therefore no rigid body acceleration is possible during the 
interval of tlffliei' ; Therefore it is evident from eqs. (2.84) 
and (2.84a) that the rigid body displacements and velocities 
at the end of the interval may be obtained by simply allowing 
the instantaneous structure to continue the rigid body 
motions existing at the beginning of the interval. Por this 
reason, it was not necessary to include the rigid body modes 
in the modal matrix in Sections 2.5.2 and 2.6.1. 

2.7 Harmonic Analysis of Big id Body Porces 

In the above sections, the force vector P in 
eq. (2.35) was treated as constant during the interval at. 
Strictly speaking, the rigid body inertia forces vaiy with 
time even within the interval at. These time dependent 
forces oan be decomposed into several harmonics with various 




amplitudes. Evidently, an improvement in the accuracy of 
the previous analysis can he achieved if the time varying 
forces are treated in their actual form specially when the 
amplitudes of the fiist few harmonics are high. Moreover, 
when the operating speed becomes near to the frequency of any 
of these harmonics, large vibrations take place due to the 
’resonance' condition corresponding to that harmonic. Such 
critical situations can not be easily detected if the force 
vector is treated as constant during' the interval. 

Unf ortimately, the harmonic analysis of the 
mechanisms is still in its infancy. As mentioned in Section 
1.1.3, harmonic analysis is available in the literature only 
for those mechanisms which have one input crank rotating at 
constant speed. Consequently, the following analysis is. . 
applicable only to that restricted class of mechanisms. 

Though the analysis may be extended to other cases, only 
planar mechanisms consisting of straight links with constant 
cross-section are considered in the following discussion. 

In this section, harmonic analysis will be done 
only for the rigid body inertia forces. The rigid body 
configuration of the instantaneous structure is considered 
as a fixed frame of reference to describe the small motion 
of the mechanism during the current interval of time A v. 
Therefore, the rotation matrix 1 in eq. (2.26), the trans- 
formation matrix f in eq. (2.38) and the matrices K, C and 
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M in eq. (2,35) are considered as time independent though, 
in a strict sense, they are also composed of sereral 
harmonics of the input angle. 

When the input crank rotates with a constant 
speed, the sine and cosine of the angles of the links can 
he expressed in the form of Fourier series. For illustrative 
purposes, such expressions for slider-crank and four-bar 
mechanisms are derived in Appendix B, Thus for any element 
E’ shown in Fig. 2.2, the sine and cosine of the angle ©g, 
can be expressed as follows; 


cosftg, 


sinSf^, 


“1 

I [cc- cos(j-l) wT + cs. sin jtoT] 
3=1 ^ ^ 


h. 


1 

[sc. cos( 3 - 1 ) wT + ss. sin jwr] 
3=1 ^ ^ 


(2.85) 


( 2 . 86 ) 


where 




to = 

T = 



total number of sine/cosine terms taken in each 
series. 

cs . , sc . , ss . = deterministic Fourier co- 

3 3 3 

efficients 

constant input speed 
time 

CS . , sc . , ss . = 0 for j > h.j . 

3 3 3 


(2.87) 



The following trigonometrical relation is used to 
find the angular velocity ug, of the element in the form of 
Fourier series. 


E * d d 

~ di~7~ ~ cos© 2 j, (sin©g,) ~ sin©g, (cos0g,) 

( 2 . 88 ) 

As the series in eqs. (2.85) and (2,86) are absolutely and 
uniformly convergent, term by term differentiation is permi- 
ssible. After substituting eqs. (2.85) and (2.86) in eq. 
(2.88), performing the differentiations and retaining terms 
up to h-j-th harmonic in the series multiplications, the 
following relation is finally arrived at. 

isi = 0 ) ^ [vc . Gos(j-1)toT+ vs. sin joiv] (2.89) 

3=1 ^ J 


wherevc . , vs. are known co-efficients. 

<) «3 

Squaring both sides of eq, (2.89) » the series 

o 

expression for Ug, is obtained in eq. (2.90) where h.^ numbers 
of cosine/sine terms only are retained. 

h. 


2 2 
tii„ . = 01 

E> 


"1 


I [nc. cos( 3 - 1 )ojT+ ns. sin jojT] 


(2.90) 


0=1 


where nc., ns . are known co-efficients. 

0 0 

The angular acceleration is found out from 
the differentiation of eq. (2.89) as follows: 



^1 

(3. 2 T* * 

®E‘ ~ dr cos(;3-1 ) UT + ts. 

3=1 

where , 

tc^ =0 

■tc^ = 3 , 3 = 2,,.., 

ts^ = -3 » 3 = 1» ••-» 1^1-1 (2.92) 

and ts- =0 hy virtue of eq. (2.87). 

The absolute normal acceleration and the 
absolute tangential acceleration of the left end M of the 
element along its axis and perpendicular to it respectively 
can be similarly expressed in Fourier series by taking 
pro3ections of the accelerations of the right end of the 
ad3acent element, (For elements generated from the same 
link, pro3ections are only formal). These accelerations 



finally take the following form; 


■^4flx “ “ ^ °o®(3-1 ) Sin j ojx] 

3=1 - 

2 V ^ 

A,, = 0) I 008(3-1 )0)^+18 . sin 3 a) 


(2.93) 

(2.94) 


where ac . , as., be. and bs . are known constants. 

3 3 3 3 

The normal acceleration and the tangential 
acceleration of the right end 1 of the element along 
its axis and perpendicular to it respectively are given by 


eqs. (2.95) and (2.96). 
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- “I. 1 (2-95) 

( 2 . 96 ) 


Upon substitution of eqs. (2.90), (2.91), (2.93) 
and ( 2 . 94 ) in eqs. ( 2 . 95 ) s-nd ( 2 . 96 ), the following series 
result : 


K 

Ux 

K 

Ny 


^1 
= a>^ I 
3=1 


^1 

1^1 


[dc. cos(j-l) (oT + ds. sin j wt] 
J J 


[ec. cos( 3 - 1 )uT+ es. sin 3 a'c] 


where , 


(2.97) 


( 2 . 98 ) 


dc . 

3 

= ac . 

3 

— 

nc . 

3 

. 1 

ds . 

= as . 
3 

— 

ns . 

3 

. 1 

ec . 

3 

=.03 

-f 

tc . 

3 

. 1 


= bs^ 

+ 

ts . 

3 

. 1 


(2.99) 


Use of eqs. (2.93) ^ ( 2 . 98 ) in eq. (2.22) yield 

the series expressions for the body forces of the element. 
On substitution of these expressions in eq. (2.25) for the 
body forces, the element oriented element load vector p in 
the form of Fourier series is obtained in eq. (2.100). 
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V 1 2 

p = [p^ cos( j-1 ) OJT + -p Sin j cdt] (2.100) 

3=1 J 

_1 _2 

where p., p. etc. are known constants. 

ti ti 

There may exist a special kind of external forces 
(i.e., the centrifugal forces of a concentrated weight, the 
inertia of a fljrwheel etc.) which have a fundamental 
frequency equal to the angular speed of the crank. These 
forces may he expressed in series of harmonics of the input 
angle and using these expressions in eq. (2.33s-) » "tiie 
element oriented element force vector p^ due to the external 
forces may be found out in the form of a series. Clearly, 
the forces like gravity etc. will contribute only to the 
first term of the series. These series expressions should 
be added to the series of p in eq. (2.100), Other t 3 rpes of 
external forces may also be present whose fundamental 
frequencies do not coincide with the speed of the crank. 

These forces, forming a separate group of series, may be 
similarly treated. However, the^e forces are not included 
in this section. 

Employing the series expression eq, (-2.100) of 
p in eq. (2.31) and noting that the rotation matrix E is 
considered as constant during the interval, the series 
expression for the system oriented element load vector ^ 
is obtained as follows. 
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p = L {pi cos( j-1 )«i}T + p^ sin jwT] 
0=1 ^ ^ 


( 2 . 101 ) 


where , 


a1 -1.-. 

p. = E p,.: ; 


^2 


T^t -2- 

R p. 


( 2 . 102 ) 


a 1 a 2 r \ 

The terms p., p., (j = 1, h. ) for all elements 

tj kJ 1 

are separately assembled for each harmonic as described in 
Section 2.4 and the system force -vector P of eq. (2.35) is 
formed as shown in eq, (2,103). 
h 

P = 


1 

1 -1 p 

I [p. cos( j-1 )wTn;+ P sin 0“'’^] 
0=1 ^ ^ 


(2.103) 


where pl and P^ are the assembled form of pi and p^ 

0 0 0 0 

respectiwely . 

As the transformation matrix T is considered 
constant during the interval, each harmonic of the input 
torque maintaining the ‘static’ equilibrium of the instan- 
taneous structure may be found out by applying eq. (2.57) 
for the corresponding harmonics of the elements of the 
system force vector P. These harmonics of the input torque 
are to be added to the corresponding elements of P in 
eq. (2.103). 

Corresponding to P in eq. (2.103), P in eq. (2.60) 
assumes the following form; 
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P _ I [ ^ j cos(3-1)w'c + 

O”*'^ 

The i-th element Pp of P may be written as 


( 


follows : 


P - = ^ L “i i 

1 3=1 


I [e.^ cos(j-1) 03*^+ fij sin 3 

i” * * * 9 


where e 


and f . . are the i-th elements of ^^P^ and $ P. 

'ij ^3 

respectively. 

When P., given by ea. (2.105), is used in eq. 


(2.63), iiie c 


omplete solution Is given by eq. (2.106). 


= e"“^' [Wj_ cos + Zi sinpiX' ] 


’ 1 + 


X [di 003(3-1)“-'+ ^3 

. A ^3 
l=i 


Consequently, 


p. (m) = 


-Opt 


[(-<TpW^ -I- |ipZp)cos lipX’ - (o^pZp + 


_ 4 


^,.W. )sar( ppxn + [<ii3 d . . 


and 11p(t) 


P. o 

^ . uifn- - 2C . “.p. 
m- 1 '3- 1 X a 


.104) 


(2.105) 


(2.106) 


sin 3 03X ] 
(2.107) 


(2.108) 
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¥/here , 


Ti? - I [d. . cos(j-1)ioi: + d^. sin ja-c] 

^ o ±2 o 


ttK? + <yt\ 


-X 

I { d^ . cos( j-1 ) a-c 
3=1 


+ dT^ sin 3 aTQ } ] 


[1 - ai e.. 

^2 13 13 0 ). 13 

i 


2(3-1) ‘■Sidi 

i__ii e. -h (1 - 


.2 2 r 

2 — SiL.) f 

1 




«ll = ° 


( 2 . 109 ) 


<3 = 


3 «i,3+1 - 


dT , =0 

i,h., 


E“i - ( 3 - 1 )“^]^ + 4 ( 3-1 


.2 2^2 ,..2 2 2^2 
(“i - 3 “ ) + 43 a a^Cj^ 


and f)^, t’ , 0 ^, Ti?, fi?, a^. and are already 


i’ t^l' '1' ' 1" 1 


def ined in eq,s . (2.63) a^id (2.69). 
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It is evident from eq. (2.109) that if the crahk 
speed CO becomes nearly equal to some sub-multiple of any of 
the natural modes of vibration, sub-harmonic oscillations 
of large amplitudes take place, especially when the damping 
term is small. Of course, sub-harmonic ’resonance' conditions 
of higher order are not important because the corresponding 
amplitudes of the force vector are usually small. Consequently, 
the amplitudes of vibrsition are not high in the presence of 
damping. 

2.8 Stress Analysis 

The dynamic equilibrium equation for each element . 
may be written in the following form; 

lEu = p + (- mu) (2.110) 

The above equation may be considered as a force-displacement 
relationship for an element where the forcing terms are kept 
at the right side. In the finite element method of analysis, 
the continuous material properties and the distributed forces 
are idealised for some selected node points and then the 
analysis for these node points are carried out as evident 
from the nature of eq. (2.35). Thus, it finally yields 
the displacements, velocities and accelerations of the 
mechanism only at the nodes of the elements, those being 



again due to the nodal forces only. If the stresses at the 
nodal points only are desired, use of eq. (2.20) is sufficient 
for that purpose. However, for any point within the element, 
the stresses are due to; (i) the actual distribution of the 
external forces acting over the length of the element, 

(ii) the nodal displacements of the element and (iii) the 
actual distribution of the vibratory inertia forces acting 
over the length of the element. lo find out the exclusive 
contributions of each one of these factors, the principle of 
superposition is used in the following manner. 

If all the nodal displacements of the mechanism 
are made equal to zero, the resulting configuration is known 
as restrained structure in the literature of structural 
analysis. The restrained structure for a slider crank 
mechanism with two divisions of the connecting rod is shown 
in Fig. 2.4. When only the external forces in their actual 
distributed condition are allowed to act on this restrained 
structure, it may be termed as 'static' restrained structure. 
The stresses Oj at any point within the element may then be 
calculated from the elementary analysis. 

Hext, all the external forces acting over the 
elements are converted into statically equivalent nodal 
forces (from the principle of virtual work) by using eqs. 
(2.25) and (2.35a). After solving eq. (2.35), the nodal 
displacements of the elements are found out as described 



77 



Fi6 S.A mirhAimo structuse 



78 


in Section 2.5.4. The stresses Oj-j within the elements due 
to these displacements alone are given hy eq. (2.20). 

At the end of each interval of time , the 
accelerations u at the element oriented element co-ordinates 
are known. Consequently, the distribution of the vibratoiy 
inertia forces may be approximately obtained from eq. (2.12b) 
as follows? 

- PAu = - PAau (2.111) 


where the left side represents the vibratory inertia forces 
at any point within the element. 

The distributed inertia forces (- pAu), as found 
out from eq. (2.111), are then allowed to act over the 
restrained structure. The stresses at any point within 

the element due to these distributed inertia forces alone 
may again be found out from the elementary analysis. Thus, 
the total stress Oj, at any point within the element is given 
by eq. (2. 112). 


^T " 




( 2 . 112 ) 


If desired, the contribution from the distributed 
element viscous forces acting on the restrained structure 
towards may be easily found out by forming an equation 
similar to eq. (2.111). The restrained structure under the 
action of the vibratory inertia forces and viscous forces 
is termed as ’dynamic' restrained structure. ■ 
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Since the shape function a is exact for beaais in 
static problems only, cr^ gives exact values. The accuracy 
of c^jjand depends on the accuracy of the modelling for 

the actual dynamic displacements, as done in eq. (2.12). 

Most of the loadings in the mechanism are 
distributed and therefore Oj and cTjjj should be included in 
the stress analysis. However, with the increase in the 
number of divisions of the links, the actual dynamic 
displacements are more accurately modelled by eq. (2.12) and 
both Cj and Ojjj decrease progressively as the restrained 
structure becomes more stiff. The displacements within the 
elements, if desired, may be obtained by the above procedure 
by using relations for the displacements in place of the 
stresses . 

In the present method, the internal stresses (and 
deflections) are calculated by completely isolating the 
contributions due to the deflections, the distribution of 
the rigid body inertia forces and the distribution of the 
vibratory inertia forces from each other. Thus, only a few 
divisions of the links and elastic modes are necessary if 
the present method of stress analysis is adopted. Of all 
available techniques, this method of the stress analysis 
requires minimum additional computations. 



2.9 Rigid Body Analysis 


2.9.1 Displacement analysis 

The kinematic analysis of a general rigid link 
mechanism may be carried out following the iterative procedure 
described in Reference (13). However, use of the transfor- 
mation matrix T developed in Section 2.5.1 may also be 
extended further for the same purpose. As the analysis of 
the rigid link mechanisms does not strictly come within the 
perview of the present study, mechanisms with re volute pairs 
only are considered here. The same analysis may be easily 
extended to mechanisms with prismatic pairs if minor 
modifications are made. 

For a rigid link mechanism with revolute joints, 
the displacement analysis implies the determination of the 
angles made by the linics with respect to a fixed line i.e., 
X-axis of the system corresponding to a specified position of 
the input links. To make the discussion simple, the same 
system co-ordinhtes are assumed at present for both the rigid 
body and kineto-elastodynamic analysis of a mechanism. An 
efficient choice of the element (and system) co-ordinates 
for the rigid body analysis is mentioned in Section 2.9. 5. 

Before explaining the procedure, it may be 
pointed out that eqs. (2.41) and (2.42) are valid for any 
combination of the lirik angles (i.e., for any distorted 
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position of the links) as long as the boundary conditions 
at the joints of the mechanism allow the rigid body motions 
(these motions are possible for any distorted shapes of the 
links) corresponding to the rigid body degrees of freedom of 
the mechanism. Thus for the improvement of some guess values 
for the link angles, eqs. (2.41) and (2.42) can not be 
utilized. Since eqs, (2.41) and (2.42) are based on the 
compatibility of the displacements j velocities etc. at the 
joints, conditions like constancy of the velocities of the 
joints in a given direction are always maintained by these 
two equations and the values of dXj^ obtained from eq. (2.58) 
for a given value of dX^ are always exact relative to the 
values of the link angles used in arriving at eq. (2.41). 

Thus, the equation for improving the guess 
values of the link angles must be nought from such conditions 
which are valid only for a rigid link mechanism (i.e., 
necessary conditions). One of such conditions is the closure 
of the loops of a mechanism. 

To describe the procedure in a very simple way, 
let a planar four-bar mechanism be taken as an example. 
Starting with some guess values and for the angles 
of the coupler and the follower corresponding to a given 
angle ©^ of the crank, the system stiffness matrix K is 
constructed and eq, (2.41) is solved for the transformation 
matrix T. For a planar four-bar mechanism, the closure , 
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equations may be written as follows; 



h 

cos©2 

4- 

cos©^ 

3 

+ 

cos©^ + 

= 0 

(2.113) 

and 

h 

sin©2 

4“ Ii-Y 
D 

sin©^ 

" h 

sin©^ 

= 0 

(2.114) 


where the subscripts 1 to 4 stand for the fixed link, crank, 
coupler and follower respectively, L denotes the length and 
0 denotes the exact angle (with the axis of the fixed link) 
of a link. The errors d0^ and d0^ in the guess values are 
given by eqs. (2.115) and (2.116). 

0^ = + d0^ (2.115) 

0^ = + d©^ (2.116) 

With the help of eq. (2.38), the errors d©^ and d©^ may be 
expressed in terms of a fictitious error d©^ (all measured 
at the level of the . guess values ©^ and ©^) in the crank 
angle ©2 as shown below. 


d9^ = 49' 

(2.117) 

40^ = 

(2.118) 


where and T„ are the elements of T corresponding to the 
rotational co-ordinates (for each link, any one of its two 
rotational co-ordinates may be used) of the coupler and 
follower respectively. Using eqs. (2.115) to (2,118) in 
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eqs. ( 2 . 113)5 ( 2 . 114 ) and solving for the following 

relations result (assuming d9^y d©^ and d©| small): 


1^ cos ©2 + 1^ cos©^ + 1^ cos©^ - 
1^ ‘Tg sin©^ + Tg sin©^ 


( 2 . 119 ) 


Ii2 sin©2 + sin©^ + 1 ^ sin©^ 
If. cos©l + 1, coioi 

j j ^ ^ 4 * ^ 


(2.120) 


Since both eqs. (2.119) and (2.120) give the 
same quantity dO^, any one of them may be used. 

As the term d©^ plays a key role in the present 
iteration process, it should be calculated as accurately 
as possible, Porthis reason, the term d©^ should be 
evaluated from eq. (2.119) or eq. (2.120) by using double 
precision arithmetic (i.e., sixteen significant digits) 
of a computer. However, it may be mentioned that the 
computations involving double precision arithmetic are 
slightly slower than calculations using single precision 
arithmetic (i.e., eight significant digits). 

The term dO^, obtained from eqs. (2.119) or 
( 2 . 120 ), may be interpreted as the error accumulated at one 
angle (here the crank angle) due to the incorrect guess 
values for the other angles. Upon back substitution of 
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this common error d9^ in eqs . ( 2 . 11 ?) and (2.118), the 
errors dO^ and d©^ are obtained , Thus , the new guess -values 
and ©)|^ for the angles may be obtained by adding the errors 
to the existing guess values as shown in eqs. ( 2 . 121 ) and 
( 2 . 122 ). 


®3 

-- Qf 

- W 3 

+ d©» 

3 

( 2 . 121 ) 

®4 

- ^4 

+ d©' 

( 2 . 122 ) 


The above process is repeated until the differences 
(or the relative difference s/between the values of d©^ and 
d©^ (or d©^) obtained from two successive iterations satisfy 
the pre-assigned convergence criterion. It may be noted that 
the iteration procedure is carried out for a fixed value of 
©2 therefore d©^ obtained from eqs. ( 2 . 119 )? ( 2 . 120 ) 
must not be added to ©2 during the iteration process. 

Equation (2.38) not only plays a key role in the 
above iteration procedure, but also provides a good estimate 
for the angles to be used as guess values .during the next 
interval of time. These estimates are obtained from the 
left side of eqs, ( 2 . 11 ?) and (2.118) when the increment in 
the crank angle (i.e., difference between the crank angles 
of the current configuration and the next instantaneous 
structure) is substituted for d©^ in the right side of these 
equations. 
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Role of eqs. (2.117), (2.118) is to transfer tlie 
errors in the gness values of the angles to the crank angle 
and redistribute again this accumulated error at the crank 
angle back to the angles of the respective links. 

2.9-2 Velocity and acceleration analysis 

Once the accurate positions of the links are 
known from the procedure described in Section 2.9.1, the 
absolute velocities at the system co-ordinates corresponding 
to the given velocities at the dependent co-ordinates are 
immediately obtained from eq. (2,38) by dividing its both 
sides by dr as shown in eq. (2,125). 

= TXg (2.123) 

• * 

where and X^ are the absolute velocities at the indep- 
endent and the dependent system co-ordinates respectively. 

Thus, the absolute translational velocities at the end points 
and the absolute angular velocity of each link are obtained 
from the appropriate elements of the vector X^. Since T 
has been obtained in the displacement analysis, the velocity 
analysis requires the trivial computations of multiplying a 
ma.trix by a vector as indicated by eq. (2.123). 

Upon differentiation of eq. (2.123) with respect 
to time, eq. (2.124) is obtained. 
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• 0 

= T Xg + T Xg (2.124) 

where X^ and ICg are the absolute accelerations at the 
independent and dependent system co-ordinates respectively. 

The term T may be obtained by differentiating eq. (2.41) 
with respect to time as follows; 

/ 

Replacement of T from eqs. (2.124) and (2.125) yields eq. 
(2.126) 

+ K^)X3 + litj (2.126) 

Since element oriented element stiffness matrix 
ic for each element is independent of the position of the 
elements with respect to the fixed frame of reference, the 

A 

matrix k for each element are obtained by differentiating 

eq. (2.29) with respect to time as is shown in eq. (2.127). 

♦ 

k = r”*^ k R -f k E (2.127) 


where E is the rotation matrix for the element. In case of 
a planar mechanism, E for an element E» (making an angle 
©T-,, v/ith the X-axis of the system co-ordinates) is obtained 
from eq. (2.26) a.s 


!i 

0 

_ 0 



(2.128) 
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where 


- sin©-p, 

COS©g, 

0 “ 


— cos9j^, 

- sin9g, 

0 1^, 

(2.129) 

0 

0 

0 



Since the absolute angular ■velocity of each link 

# 

is known from eq. (2.123), the matrices k for each element 

is known from eqs. (2.23), (2.12?) and (2.128). Assembly of 

♦ 

the matrices S; for each element may be done by using the 

integer ma.trix NS as described in Section 2.4 to form the 
• • * 
system matrix K and consequently the submatrices and 

Using the matrices and the given values of X^, X^ 

in eq. (2.126), the vector X^^ is obtained. The absolute 
angular acceleration and the absolute tia-nslational accel- 
erations at the end points of each element are known from 

• • 

the appropriate elements of the vector X^. The absolute 
velocities and accelerations of each element in its axial 
and transverse directions may then be easily found out since 
the angles of the links are also known. 


2,9.3 Quasi-static analysis 

To apply the method of analysis described below, 
it is necessary to define separate system co-ordinates 
(as shown in Uig. 2,3a and Pig. 2.3b) along the directions 
of the bearing reactions with numbers following that of 
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the last dependent system co-ordinate. Let the combined 
group of the dependent co-ordinates and the co-ordinates 
corresponding to the bearing reactions be denoted by the 
subscript 0. The group of co-ordinates corresponding to the 
independent co-ordinates is denoted by the subscript A as 
was done in Section 2.5,1. 

The instantaneous structure of a rigid link 
mechanism is in equilibrium under the actions of the external 
forces, the rigid body inertia forces, the variable input 
forces (calculated from eq. (2.57)) and the bearing . reactions . 
Since each co-ordinate corresponding to the bearing reactions 
represent an additional rigid body degree of freedom (which 
is suppressed, in actual case, due to zero displacements at 
these co-ordinates), the force equilibrium equation for the 
rigid link mechanism given by eq. (2.57) may be extended to 
find the bearing reactions as follows: 

Tq + + T^ = 0 (2.130) 

where the vector Tq contains the variable input forces and 
the bearing reactions, and are the rigid body forces 
(external forces and the rigid body inertia forces) at the 
group of co-ordinates C and A respectively. The matrix T^ 
in eq. (2.130) is found out from eq. (2.131) which is an 


extension of eq. (2.41). 
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T 


K” TC 
■ ^AA ^AC 


(2,131) 


It majr' be noted that the first set of the columns of 
contains the matrioc T defined by eq. (2.41) and the first 
set of the elements of the vector contains the variable 
input forces. Thus, elimination of T^ from eqs. (2.130) and 
( 2 . 131 ) yields the variable input forces and the bearing 
reactions, contained in the vector T^, as shown in eq. (2.132). 


Tq - - E'q + 


( 2 . 132 ) 


-1 

It may be noted that the static deflections of the 

mechanism due to the rigid body forces are obtained from 
eq. ( 2 . 132 ) as an intermediate product. Equation (2.132) 
may also be derived from the force-displacement relationship 
of the mechanism. 

Once the bearing reactions are obtained from eq. 
( 2 . 132 ), the rigid body pin forces at the moving ends of the 
input links may be obtained from the equilibrium of forces 
acting on the input links. Erom the compatibility of the 
internal forces acting at these moving ends, the pin forces 
at the same ends of the next links connected with these 
joints are found out. In this way, the pin forces acting 
at all the joints of the mechanism may be easily found out 
by successive applications of the equilibrium equationsfor 
the links and the compatibility equations for the joints. 
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When more than two links are joined in several pin-connections, 
the number of simultaneous equations to be solved may be 
considerably reduced if the above process is initiated from 
all support-ends at a time and if provisions are kept ?/ithin 
the computer program to choose the successive links and 
joints judiciously. 

As mentioned in the beginning of this section, 
it is not necessary to use the same system co-ordinates for 
the rigid body analysis and the kineto-e las to dynamic analysis. 
For example, divisions of the links are not necessary for 
the rigid body analysis and also for determining the quasi- 
static deflections. Moreover, for a link with hinged joints 
at both ends (i.e., a truss), only translational co-ordinates 
(two co-ordinates at each joint for a planar link) are needed 
to describe its rigid body as well as elastic behaviour. 
Therefore these types of co-ordinates should be used for the 
rigid body analysis and quasi-static deflections of the 
mechanism to minimize computations. The shape function and 
the generation of the rigid body forces for these type of 
co-ordinates are given in Appendix D. This shape function 
is exact for the quasi-static analysis; therefore the nodal 
deflections obtained from eq. (2'. 132) are also exact. The 
quasi -static deflections within the links may be calculated 
following the procedure of Section 2.8 which are exact for 
an exact shape function. One difficulty in using this type 
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of co-ordinates lies in the calculation of the errors in the 
angles (eqs. (2.11?) to (2*120)), to be calculated 
indirectly from the errors in the translational co-ordinates. 
To circumvent this difficulty, one additional co-ordinate 
may be used for rotation of each link and’ the associated shape 
function may be easily determined following the procedure 
described in Section 2.3. 

Por a planar four-bar mechanism, the rigid body 

displacements, velocities, accelerations, pin forces, variable 

input torque and quasi-static deflections have been found out 

(using shape function of eq. (2.13)) for a complete cycle 

by means of the present procedure with an increment of 

4 degrees in the crank angle and a convergence criterion of 
-3 

10 . The entire analysis takes only 0.1 min. in IBM 7044 

and the iteration process converges in three iterations (on 
the average) using single precision arithmetic. This computer 
time is much less compared to the existing matrix analysis 
procedure^ 

Like the existing matrix procedure, the analysis 
procedure described above is completely general and is ideally 
suited for a digital computer. In addition, its main advan- 
tages are: (i) nximber of computations is less by an order 
of magnitude, (ii) gives quasi-static deflections (and 
therefore stresses) as an intermediate product, (iii) nay 
be applied to mechanisms with multi-degrees of freedom, and 
(iv) simple dimensions of the curved links are needed. 



CHAPTER III 


EYHAMICS OE AXIALLY LOADED MOVING LINKS 

3.1 Introduction 

When a mechanism operates, all links of the 
mechanism experience varying translational as well as angular 
accelerations. Moreover, the axial accelerations and the 
interactions "between the links transmitted through the joints 
generate space dependent and time dependent axial forces 
within the links. At low operating speeds, the effects of 
these accelerations on the vibrations of the links are not 
significant. However, even in this case the axial forces 
transmitted through the pins and their effect on the d 3 mamics 
of the links may be considerable if the magnitudes of the 
external forces acting on the mechanism are high. Since the 
accelerations on the links increase approximately with the 
square of the input speed, their role on the motion of the 
links becomes increasingly significant at higher speeds. In 
the previous analysis described in Chapter II, these 
accelerations and the external forces were considered to 
contribute towards the mechanism behaviour exclusively in 
the form ef rigid body forces. Once the force vector was 
formed, these accelerations were considered to be temporarily 



absent and the vibration of the instantaneous structure due 
to. this force vector only was analysed during each interval. 
In addition, the effect of the axial forces of any kind on 
the tramsverse vibration of the links was not included. 

The axial forces- acting on a link may ariss due to 
the following reasons; (i ) rigid body axial accelerations 
within the link, (ii) the pin forces transmitted by the other 
links through the joints, (iii) the elastic axial acceler- 
ations within the link, and (iv) the elastic forces generated 
from its axial shortening. The axial shortening of a link 
may be due to the axial forces or the transverse forces or 
both acting on it. In this chapter, the accelerations of 
the links are considered to be active during the interval 
of time At and the vibrations of the links in these acceler- 
ation fields are analysed. In addition to this, the effect 
of the axial forces on the transverse vibrations of the 
links, which is essentially non-linear in character, is also 
included in an otherwise linear analysis. To achieve this 
objective, the equation of motion is derived in Section 3.2 
for a general case by first using a rotating co-ordinate 
system and then Lagrange’s equation of motion. The inclusion 
of the effect of the acceleration field and the axial forces 
on the vibrations of the links results in an antissnnmetric 
secondary mass matrix and a symmetric geometric stiffness 



matrix respectively. The solution, of the equation of motion, 
with these antisymmetric matrices as coefficients, by the 
normal mode method requires the use of the complex eigen- 
values and eigenvectors. To avoid the complex arithmetic 
and the accompanying excessive computations, an approximate 
method is outlined in Section 3.3 which solves the equation 
of motion with sufficient accuracy but using lesser amount 
of computations. When the elastic axial forces are considered, 
the geometric stiffness matrix becomes non-linear, lor 
this case an iterative procedure is presented in Section 3.4, 
at the end of this chapter, 

3.2 Derivation of the Equation of Motion 

In this section, the equation of motion will be 
derived for an element E» in a more general situation, i.e., 
when the element is axially loaded and is subjected to the 
longitudinal and transverse vibrations in an acceleration 
field as shown in Fig. 3.1 . Choosing a fixed co-ordinate 
system X-Y, let the co-ordinates of the left end M of the 
element E’ at any instant be X, Y and its inclination with 
the X-axis be Q^, as shown in Fig. 3.2. 

At this instant, the absolute rigid body 
angular velocity and angular acceleration are ojg, and Og, 
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respectively and the axial forces at the left end M and 
right end F are and P^ respectively. A rotating co-ordinate 
system x-y is chosen along the axial of the element and 
perpendicular to it which remain fixed with the element and 
rotates along with the element, e ^ , "Ig and e^, are the 

unit vectors in the x-y and X-Y co-ordinate systems respec- 
tively (Pig. 3.2) satisfying the following relations; 



= 


- 

-4“ 

II “ 

(5.1) 

® "1 • ®2 

= 

^2*^1 



= 0 

(3.2) 

CD 4^ 

H 

11 


c osQgi f 

-4- 

- eg 

sinOg, 


(3.3) 

- 

-4- 

: e^ 

sin©g 

1 + ^2 

COS&gj, 


(3.4) 

» 

e-j = 

“b’ 

p 0 
^2 ^ 

•4- 

eg = 

- 


(3.5) 

= 

->■ 

-4* 

= e 

I 

011 = ® 


(3.6) 

a 

X = 

X = 

• 

= y = 

^ 7 = 

0 


(3.7) 

Let 

any 

point 

within 

, the element. 

, initially at 

the 


position J with co-ordinates x, y, take the position K with 
co-ordinates x + u^, y + u^. (all measiu-ed in x-y system) 
due to the longitudinal and transverse vibrations of the 
link. The final position vector Eg- is given by the following 
equation; 
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Eg. = X ej + Y Bj-j. + (x + u^) + (y + u^) (3.8) 

Differentiating eq. (3.8) with respect to time and using the 
relations given by eqs. (3.3) to (3.7), eq. (3-9) is obtained. 

^ + 1 sinOgj + — wjji (y *^1 

• ® • »i 

- X sinSg, + Y coseg, + u^ + oig, (x + u^) (3.9) 

The kinetic energy of the element is; 

• * 

Kg = ^ J p Rg.:^ dv (3.10) 

where p is the mass density of the link and the integration 
is carried over the whole volume v of the element. 

Introducing eq. (3.9) in eq. (3.10) and utilising 
eqs. (3.1) and (3.2), the following relation results. 

Kg = ^ ! P [(K cosOg, + Y sineg, - <<^(7)^ + (* + 

• • 2 * 2 
+ (- X sinBg, + Y cos©g, + u^.x) + (u^ - £i^,Uy) 

® • m 

+ 2(X cosGg, + Y sin0g, - «£ty)(n^ - “s’^y^ 

+ 2(- X sinOg, + Y cosGg, + ajg,x)(Uy + ajg,u^)]dv 

(3.11) 

The continuous displacements n^(x, y,v) and 
u (x, y, t) may be replaced by the nodal displacements u..^ to 

y 

Ug (Rig. 3.1) in the element oriented element co-ordinates 
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x-y with the aid of the shape function a defined by a 2x6 
matrix in eq. (2.13) and is rewritten in expanded form as 


follows: 


.1 


3 ^3 


= .L ^23 "^3 


(3.12' 


( 3 . 15 ) 


Employing eqs. U* ‘ 
"b63?ms S.P6 fomid oni 


. (3.12) and (3.13) in ©I- (3.11), the following 


1 = / P [(- i sin©g, + Y cos©g, + '*fetx)«E'^ 1 j 

I , y 

^ _ (i cos%,, + i sinGg, - t^,y)'^’^23 


K - “E'^y^“E '^23 




= / P [(X cos©g, + Y sine^t - 

5u. V 

3 

+ (- X sinSg, + Y coseg, + 


The following relation is derived from eqs 
and ( 3 . 15 ) for further use. 


. ( 3 . 7 ), ( 3 . 14 ) 



3K 

d /__e 


( -“ ) - -T 




hlx * - 


^ LV + ‘"E 








2o3-p I li.^ 

+ /■ P [^pj ~^l i ^ ‘ 

V 23 13 |_2Bj,,UyJ 


/p [Spj -a^j] 


a-ptW-y 
® ^ dv 


L“b«\ 


(3.16) 


where » 


A, - X cosa-P, + f Sin%, = Acceleration of the point M xn 

A ^Yie x-direction 


X sin« + Y cosep, = Acceleration _ of the point 
-X sin«g, ^ M in the y-dxrection 


= ©TPt 

E» E’ 

«• 

Et = ^ ®B’ 

The potential energy stored in the link E« 


(5.18) 

(3.19) 


( 5 . 20 ) 


(neglecting the gravity) is the 
is given hy eq. (3»21). 


e strain energy of the link and 


pi; e dv 


(3.21) 


where 


Ml 


a(xy y, "c) and e(x, y, t) ere the stress and strain 

at any paint of the link. 

Continuous functions c and e can be expressed in 
terms of u as shown in eqs. (2.20) and (2.16) and are rewritten 
hers, in- tiie following expanded form. 


e = 


(3.22) 


where , 


I. 




ill ’’23 "3 

ti 


(3.23) 


e 

A ®2 


(3.24) 


(b, E, G, A and A are 


defined in Section 2.3) 


From eqs. (3.21) and (3.22) 


/ (a^e^ + 


(3.25) 


Using eqs. (3.23) and 

P with respect to u-j 

6 


(3*24) (3*25) Q-tid dif f sx' 611 'fc 


= 5 i - V23 " — hj ^2)^" 

GA 

= / (B .e, + -J- '>23*23'^'^ 


I *^23 ^ 


1 0 
GA. 
0 -X 


(3.26) 



The strain energy given by eq. (3.21), when used 
along with equation(2.17) for strains, does not take into 
account the additional work done by the axial forces on the 
foreshortening of the link due to the transverse displacements 
i.e., the buckling effect of the axial forces. Por the time 
being, it is assumed that these axial forces come from the 
rigid body forces (i.e., the rigid body axial accelerations 
of the link under consideration and the rigid body pin forces 
exerted by the other links on it) alone and therefore have 
known values. The additional axial forces generated due to 
the elastic deformations of the links are considered later 
in Section 3.4. 

To include the buckling effect in the finite element 
method of analysis, the conventional procedure^ ^ is to use 
the higher order terms in the expression for the strains in 
eq. (3.21). But this procedure is unable to take into 
account the variation of the rigid body forces along the 
link, For this reason a special work function^ ^ is intro- 
duced below. 

The displacement al° of the point of application 
of the axial force acting on the element (+ve in the +ve 
X -direction) due to the transverse displacement Uy alone is 
approximately given by eq. (3.27). 



(3.27) 


>, . Su ^ 2 

% 2 


Denoting the rigid body axial force at any point x by 1 , 

I5C 

the work done by the foreshortening of the element due to 
the bending is given by eq. (3.28) 


where , 


and 



X 


W. 


1 

2 


/ P, 


9u 2 
(-vr^) dx 


X ^ 3X 


(3.28) 


the work done 

the rigid body axial feme acting at the left 
end of the element dx (considered +ve in the 
+ve x-direction) . 


The axial force P , shown in Pig. 3.3s is* 

IX. 

X 

^x ~ (3.29) 

O 

The rigid body pin force may be obtained either by direct 

analysis of the rigid link mechanism (as described in 

Appendix D) or by the general method presented in Section 2.9~. 

3U 

Prom eq. (3.13), oa.n be written as 


3X 



(3.30) 


e ' 

®' 2;3 


9a 


22 . 


9X 


where 


3 = 1 > • • • > 


6 


(3-31) 






Introducing eq. (3.3©) in eq. (3.28) for — ^ , eq, (3.52) is 
obtained . 




an. 

J 


= / 


O 


3U 

I a ’ . ■— 
X 23 8x 


(3.32) 


The Lagrange's equations of motion for the link 
can be written in the following form: 


dt 


aK^ 9K^ 

as. ail. 

J J 


ap 


au . 


aw, 


an. 

J 


= 


( 3 . 33 ) 


j = 1 , 


where . is the external force (if any) acting at the j-th 

t] 

element co-ordinate. 

Writing eqs. (3.16), (3.26) and ( 3 . 32 ) for j = 

I 5 2 , ... 6 , introducing them into the set of equations given 
by eq. (3.33) and using eqs. ( 3 . 12 ), (3.13), (3.23), (3.24) 
and ( 3 . 29 ), the following equation of motion in matrix form 
is finally arrived at. 


/pa' 


V 


L^y + “E'^ “b’^J 


dv + m u 


+ O',,-,, m^ u + E u + u = Q 


£j 

where , 


m = / pa^ a dv 

V 


2 - 

a>g,m 


u + 2 m 


a 


“E' 


u 


(5.34) 

(2.24) 



liiG 


[m^.] 


(3.35) 

E = 

J b^ X 'B dv 

(2.23) 


Y 



= -| ^2± ^2) 

(3.36) 

■0, 

= ■f , . . . , } 

(3.37) 


The matrices 


a , mG 

m and K 


in final form are given in Appendix B. 


If the x-axis passes through the centroid of the 
cross-section of a straight link, the following relation 
holds . 


/ y dv = 0 


V 

Noting that the second row of the shape function does not 
contain y, the first term becomes equal to 

_ -a 



(3.38) 


where , 



= -/pa' 

T 


■“Eiy 


dv 


(3.39) 

(3.40) 


and 


_a 

P 


o 



I>i7 

It may be noted that p in eq. (3.39) is same as 
that given by eq. (2.25) differing only in notations and 
p®' is the additional rigid body force due to the effect of 
the angular acceleration Og,. This term is small when either 
ttg, or the second moment of area of the cross-section is 
small. 

Thus the general equation of motion (3.34) is 
rewritten as follovre; 

ffi !i + 26^,# d + (k - a)|,m + k^ + aglffi®)u = p + p® + ^ 

(3.41) 

Mi 

The terms m u, k u and p were also previously 
obtained from the conventional procedure used in Chapter II. 

cX * 

The term 2a3g,ffi u comes due to the Coriolis component of the 

elastic longitudinal and transverse acceleration of the link 

2 

in a rotational velocity field. The term -(i>g,m u takes into 
account the centrifugal effect of the elastic motions (in 
both directions ) in a rotational velocity field. The 
additional acceleration term arising due to the longitudinal 
and transverse vibrations in a rotational acceleration field 
is given by ag,m®u. The term u includes the buckling 
effect of the rigid body forces as explained earlier. It 

' ' 3 » 

is interesting to note that the secondary mass matrix m , 

given by eq. (3.35), is ant iS 3 nmne trie while the other matrices 

■ ; O' ' ' ' ' 

m, k and E are symmetric. 



For the reasons explained in the next sections, 
the element matrices and the force vector in eq. (3.41) are 
divided into five groups, i.e., (i) m, (ii) 2aig,m®', (iii) 

(E - <jjg,m + k ), (iv) agjSl and (v) the vector p^ = p + p^ + Q 
The first four groups of the matrices for all elements are 
transformed into the system oriented element co-ordinates and 
next are assembled to fom four system matrices, following 
the procedure described in Sections 2.3 and 2.4* It may be 
noted that the pin forces % are the internal forces and 
therefore need not be included in the force vector because the 
stiffness (or displacement) method of finite element analysis 
automatically satisfies the compatibility of the internal 
forces whenever the compatibility of ttB displacements are 
satisfied. 

When the links of the mechanism are subdivided 
into several elements, equation of motion for each element 
is evidently given by eq. (3.41) where the matrices of eq. 
(3.41) are to be formed for each element. For a link with 

3 , O 

constant cross-section, the terms ffi, 2ii)g, SI. , k, gt ni? 
a^, m and p are identical for all the elements generated 

from that link and are to be calculated only once . However, 

n — 

the terms k , p and Q are different even for these elements 
and therefore have to be calculated separately. 
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The assembly of the element matrices and the 
incorporation of the special conditions of the mechanism into 
the system matrices are carried out as described in Section 
2.4. The general equation of motion in the system co-ordinates 
finally assumes the following fom: 

m + (C + C^)U + (K + + K^) + U = P (5^42) 

where, 

M = mass matrix 

C = damping matrix 

C = apparent damping matrix due to the additional 
Coriolis components of acceleratior^ 

K = elastic stiffness matrix 

= apparent stiffness matrix due to the additional 
normal components of accelerations 

Q. 

K = geometric stiffness matrix 

T 

K = apparent stiffness matrix due to the additional 
tangential components of accelerations 
P = force vector 

and UCt) = displacement vector in the system co-ordinates. 

The matrices M, C^, K, K^, and are formed 

£l 

from the assembly of the element matrices m, 2a3g,m , k, 

and of eq. (5.41). Similarly,, the force 

vector P results from the assembly of the element force 
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vectors of eg., (3.41-). The matrix C is the damping matrix 
of the system and may be taken (as was done in Section 2.5.4) 
to be proportional to the critical damping at the, normal modes 
of the mechanism. Though associated with the damping term, 

Q* 

the matrix C is not dissipative because the contribution of 
this term towards the total amount of work done by the 
individual normal modes of the system is zero. The conser- 
vative nature of 0 follows from its antisymmetry. 

3.3 Solution of Tkjuation of Motion 

The undamped eigenvalue problem for equation (3.42) 
may be formulated in the form of eg. (3.43). 

MX + + K°X = 0 (3.43) 

The eigenvectors (or eigenmodes) obtained from the solution 
of eg. (3.43) will form a system of uncoupled equations in 
the normal co-ordinates. But these eigenvectors and the 
corresponding eigenvalues are complex quantities and excessive 
computations are involved with their use. 

Alternatively, the eigenvalue problem may be 
formulated in the following way; 

MX + K°X =0 (3.44a) 

where = K + + K (3.44b) 



Ill 


Because of the antiS 3 nDiiietuy of K , the eigenvectors 
of eq, (3,44a) are mutually independent but not orthogonal to 
each other with respect to either M or Pox this reason, 

the solution of eq. (3.44a) by Power method to find the 
sub-dominant eigenvalues and vectors requires ^ either the 
right-hand and left-hand eigenvectors simultaneously or the 
matrix deflation procedure. Either of these two techniqi;!es 
increases the volume of computation. Moreover, the eigen- 
vectors obtained from eq. (3.44a)can,not diagonalise the 

Si 

matrix C and so the resulting equations in the normal 
co-ordinates still remain coupled. 

In view of the above, the following simple formula- 
tion of the eigenvalue problem is adopted for subsequent use. 

m + A = 0 (3.45a) 

where = K + + K^, the symmetric part of (3.45b) 

As both M and are real and symmetric, the eigenvalues and 
eigenvectors of eq. (3.45a) are real and orthogonal to each 
other with respect to M or K®. But these eigenvectors do not 
diagonalise K®" or C^. lo overcome this difficulty, an 
approximate method is developed below for solving the coupled 
equations in the normal co-ordinates. 

Two assumptions are made in the following 
approxinate method. The first of them say be stated as 




follows; the eigenvalues and eigenmodes obtained from eq. 

( 3 . 45 ) represent the actual behaviour of the system (in its 
free vibration) with reasonable degree of accuracy. The 
second assumption will be mentioned in due course. Though 
the approximate method described below is equally applicable 
to all situations treated in Chapter II, it is assumed, to 
simplify the presentation, that the rigid body forces remain 
constant during the interval of time at and satisfy the 
static equilibrium equation for the rigid body degrees of 
freedom (so that only the elastic modes may be used). 

The eigenvalues and eigenvectors of eq. (5-45) are 
found out following the procedure described in Section 2.5.2 
and the modal matrix is formed. Employing this modal 
matrix $' , the following transformation of co-ordinates is 
defined. 

U = (3.46a) 

where T)(t) is the displacement vector in the normal co-ordinate 
system. As is independent of time (during the current 
interval), the following equations result from the differen- 
tiation of eq. (3.46a) with respect to time t. 



C = 

(3.46b) 

and 

11 = 

(3.46c) 
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Substitution of eqts. (3,46) in eq. (3-42) and pre-multipli- 
"ti 

cation by $’ , yield the following equation of motion in the 
normal co-ordinates. 

Mtl- + (C + C^)in + (K® + = P (3.47) 

where , 


1 = M 

C = G $' 

(3.48) 

and P = P 


g 

As explained in Chapter II, M, C and K are 
diagonal matrices. Consequently, the i-th equation of eq. 
(3.47) may be written as; 

m m 


'll + 2li»i>1i + MiT,. + = 


P. 


(3.49) 


3=.i. 


r — ■ '1 , . . . , 2]^ • 


where p^ is the i-th element of P, cu^ is the i-th natural 
frequency found out from eq. (3.45), nj_ is the i-th element 
of P, 5- is the damping ratio at the i-th mode, and v^- 

are the (i, j)-th element of C®" and M~^K^ respectively* 
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.a 


Since C , K are antisymmetric, 


11 


'ii 


= 0 


(3.50) 


Consequently, it is convenient to rewrite eq. (3.49) in a 
rearranged form as shown in eq. (3.51) 


h + “1% = Pi - 


m 

y e ^ .T] . 

ail ^ 

a^i 


m 

A 

aVi 


"a 


(3.51) 


i = 1 , . . . , m. 


In view of eq. (3.50), the left side of eq. (3.51) still 
maintains the standard form of eq. (2.63) and shows that the 
coupling terms affect the forcing terms more significantly 
than the natural frequency oj^^. Thus, the first assimption 
statei above can be further reinforced by adding very few 
additional modes in 

The complete solution of eq. (3.51) can be divided 
into three parts: i) the homogeneous solution ii) the 

particular solution pT due to the forcing term p^ only and 


iii) particular solution P^ due to the additional terms 


m 


given by the following relation: 


Thus, the complete solution is 


p^(v) = p^('*^) + 'n?('^) + ^j_(^) 


T < T < i: + A T 

o o 


(3.52) 
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where , 


r\. 


- 


1 2 

cos{ij^T’ + sinii^iE' ) 




m.(o. 


-f — 


% - 


(3.53) 


a . 
1 


P'i 


(1 - ?2)Vi 


^ = damping ratio at the i-th noimal mode 
1 2 

I ,1. = constants of integration 

and is given by the particular solution of the following 

equation; 


m 


ili + + „|n| = - 


(3.54) 


(the restriction 3 i is removed in view of eq. (3.50)). 

Since the right side of eq. (5.54) contains various 
functions of time x, the functional nature of the actual t)^ 
is vezy much involved. Therefore, instead of finding the 
actual (which satisfies eq. (3.54) identically), an 
approximate expression for is sought which satisfies eq. 
(3.54) only at some specified instants of time. The second 

Q. 

assumption is now made as follows; the actual satisfying 
ei. (3.54) can be adequately described by a parabolic 
expression (given by eq, (3.55)) during a small interval 


of time AT. 




r)l{x) 


3 1-2 3 !■ 


•" -^o)' 


(3.55) 


where the unknown constants ^3i-l ^31 

determined. 

As 'c = X- is small and t)^ is expected to 

vary widely in this small interval ax, the second assumption 
S66I11B "bo bo jus'tifiBd, 

Eq. (3.52) can now he rewritten in eq. (3-56) as 


. (t;) = llfbx) + ^ ^3i-2 ^3i 

^ d .1 m /.v. 


.2^2 


r] 


Hm- 2 + -■'o) 


m. 0). 

11 


+ Hj.Ct - 


(3.56) 


consequently , 


hj_(x) 


llf?(x) + + ^3i-1 ^31^”^ 

11 J- ^ 




(3.57) 


where 


f;('^) 

f^(x) 


fi(v) 


-CiX’ 

e cos pj_x' 


- a . x’ 

e ^ sin p.. 


(f]_(t)) = e ■’i t'r' 


aiT 


s in IJ. 'c * ) 

(3-58) 


P-i '"i 


4(-t) 


0^- X 


1- (f^(v)) = e ^ (- aj_ sin 

i 

+ p. cos Pj_x' ) 


and 
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Using the initial conditions 


tlifr^o) 


r\. 


Ui ( Xo ) = 

-1 ^2 


(3.59) 


the constants of integration are found out from eqs 

(3.56) and (3.57) as 


4 


4 


(hi - ' 2 “ ^3i-2^ 

m-o.. 

ri? - ll" + a.i] 
ji 1 1 1 . 


(3.60) 


Eliminating and I? from eqs, (3.56) and (3.57) 
by eq. (3.60), Tij_(T) and Ti^(x) are reduced to 


rii(T:) = pj_(T) + %i.-2 ■^75i-'l(x) + Y3i(x) 

(3.61) 

fiij_(T) = 9j_(x) + ’l^3i_2('^) ^3i_2 ’^31-1^'^^ ^3i-1 %i 

(3.62) 

where , 


Pi ^.1 


P^t) = {r,° - + [^i + -Ji (q 


m. 0) 
11 


r)3 


f.(x) P, 




Y3i_2('c) = 1 - 4(x) -^fi(x) 


2^- Pi 






31-2 
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Y 


3i-1 


('c) 


T - 

0 




Yj.(t) = (x - x^)' 


9j_('c) = (t 1? - — ^)q('i:) + [r)° + a^(fi° 


m.o) . 
1 X 


p. 


m-uY" t*! 
11 


(t) = - fj(i;) - 


a,ft(x) 


1 1 


it(t) 


♦31-1 1-') = 1 


(3.63) 


4-3^ (T:) = 2 {x - X^) 


Substituting ) , f)^(x) and rij(x) fromeqs. (3.61), (3.62) 

and (3.55) in eq. (3.54), the following equation is obtained, 


4 ^ 51-2 + - •'o + + [(x -x„)2 + 4Ci»3_(x - x^)] H 


m 

= - . 1 ^ [^iD^ ^ '<'30-1^'^ ^^33-1 

+Y5^(x)H^^>3 (3.64) 

Since eq. (3.55) contains three unknowns, three 
values of X are to be used on both sides of eq. (3.54), i.e. , 
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eq. (3.64) to determine these unknowns. But ^ eq. (3.64) 
is coupled and hence all the unknowns 3 = . . . , 3 m are 

to he found out from a set (3m in number) of simultaneous 
equations . 

The following three values of f are used 




= t: + 
0 


A'C 


(3.65) 


1 2 

Alternative choice of t°, r and % could as well 
be made by using Chebyshev's polynomials-^^^ ^ Bor each value 
^ ...j m, the above three values of are substi- 

tuted in eq. (3.64) resulting in 3 simultaneous equations. 
Finally, the 3m noimbers of equations, thus obtained, are 
rewritten in the following matrix form. 


&H = L 


( 3 . 66 ) 


where , 


^3i-q 3=1 " '^3 


1 '’ij 


® 31 -q. 33 -r = Uj 'P33-r(-'‘'^ ^ '’ij P^33-r(^’’ 


(3.67) 


3 ^ i 


®3i-q,3i-2 


2 

“i 


Si-q,3i-1 

^3i-q,3i 





i 

3 


1 , - . . , m 

j ...» ^ 
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q is given values 0, 1 and 2 successively for each, value of i. 
Similarly r is given values 0, 1 and '2 successively for each 
value of j . 

Since L and G- in eq. (3.66) are known, the unknowns 
H are obtained by inverting the matrix G as follows; 


H = G“''l 


( 3 . 68 ) 


Using these values of H in eqs. (3.61) and (3-62), 
the displacements p and velocities p in normal co-ordinates • 
are obtained. The accelerations 11 in normal co-ordinates are 
obtained from eq. (3.51) as follows; 


Hi = Pj 


m 


'13 


f). + V. . p.) - 2C. T). 

U 13 0 1 1 '1 


2 

Th . 
1 1 


(3-69) 


The displacements, velocities and accelerations 
are next found out in the element oriented element co-ordinates 
following the same procedure as explained in Chapter II. 


3.4 Elastic Axial Forces 


In the above sections, the pin forces in eq. 
(3.36) were considered to be independent of the elastic 
behaviour of the mechanism and were evaluated from the 
analysis of the rigid link mechanism. In actual case, 
additional axial forces within a link are generated due to 
its own longitudinal vibrations, the foreshortening due to 
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its transverse vibrations and the elastic effects of other 
links transmitted through the pins at its ends during the 
interval at. In other words, the term S' in eq. (3.36) is 
a function of the element displacements (represented by ^(t)) 
and therefore eq. ( 3 . 42 ) becomes non-linear. As a good 
estimate of the axial forces is available from the analysis 
of the previous interval, the iterative procedure, as 
described below, seems to be more suitable than the other 

methods available in non-linear finite element method of 

, .. (58,60,65) 

analysis^ ’ ’ . 

At the beginning of the current interval, the 
element oriented element displacements (axial) u^(t^) and 
u^(Xo) at the two ends of each element are known from the 
analysis of the previous interval. Then the average value 
of the axial force E acting on the element at that instant 
is found out from eq. (3.70). 




(3.70) 


where E, A and 1 are Yoimg's modulus, cross-sectional area 

and length of the element. Evaluating!^ in this way, the 

Or 

element geometric stiffness matrix E is formed from eqs. 
( 3 . 29 ), ( 3 . 31 ) and ( 3 . 36 ). These matrices are next assembled 
to form the geometric stiffness matrix for the whole 
mechanism and is ixsed in eq. (3.42), After solving eq. (3.42), 




oo 


the element displacements (t + Ax ) and u, (t + At) at 
the end of the interval is determined as described in the 
previous section. Using these element displacements, the 
axial forces for the elements are calculated from eq. (3.70). 
The difference iE betv/een this value of F and its foimer 

-A. UC 

value is evaluated for each element. If the maximum of 
these differences is less than the prescribed value, the 
procedures for the next interval should be started. Other- 
v;ise, a new geometric matrix for each element is found by- 
using the latest F in eq. (3.3i) and all such matrices are 
assembled to form a new geometric matrix for the whole system. 
Employing this new geometric stiffness matrix, eq, (3.42) ' 
is again solved for the element displacements and eventually 
for axial forces. The process is repeated for the same 
crank position until the convergence is achieved. The 
element displacements, velocities and accelerations at the 
end of the iteration process are used for the stress 
calculations of the current interval and as the 'initial 
conditions' for the next inter-val as mentioned in Chapter II. 
Since at the start of the iteration, elastic axial forces 
are found out from eq. (3.70) by using the values of the 
displacements obtained from the previous interval, the 
calculation for the rigid body axial forces may be omitted 


altogether. 



123 


Evidently, the above procedure is more time 
consuming than the equivalent procedure for the rigid body 
axial forces because eq, (3.42) (which requires the solution 
of the eigenvalue problem, eq. (3.45), and eq. (3.66) in 
particular) is to be solved for each iteration. However, if 
the interval of time is made smaller, very few (usually t?/o 
only) iterations are needed for each interval. 

It may be noted that though the matrices K and K 

s 

in eq. (3.45) are singular, the matrix K is, strictly 
speaking, non-singular due to the presence of the non- 
singular matrix K^. However, the elements of are small 
compared to those of K and K (except at unusually high 
speeds). Thus, the matrix K is always nearly singular. Eor 
this reason, it is advisable to use the method of Section 
2.5.2 for the solution of the eigenvalue problem. 



CHAPTER IV 


DIRECT S!EEADY STATE AMLISIS 
4.1 Introduction 

The basic framework of the analysis presented in 
the Chapters II and III is based on the conventional procedures 
of the finite element method. Though this framework allows 
the designer to find the deflections and the stresses of the 
mechanism in very general cases, it nevertheless fails to 
predict the parametric instability of the mechanisms unless 
frequency response curves are obtained at the cost of huge 
computq? time . Moreover, the method of analysis treats the 
continuous motion of the mechanism by dividing it into 
numerous small intervals, the vibrations at the end of each 
interval being incorporated in the subsequent interval as 
’initial conditions'. Thus, even if steady-state deflections 
only are desired, the analysis has to be carried out 
through a large number of intervals until the steady-state 
conditions are reached. This method of describing the 
steady-state motions by recurrent use of the transient motions 
is certainly not efficient as it entails a heavy loss of 
computer time. In view of these drawbacks a new approach 
to the analysis of flexible mechanisms is presented in this 



chapter. This method of analysis, applicable to a restricted 
class of mechanisms for which the rigid body harmonic analysis 
is available beforehand, enables the designer to find out 
the steady-state deflections and stresses in a straight- 
forward manner. The analysis may be equally applied to 
determine the instability zones of the ope 2 rating speeds of the 
input crank. In Section 4.2, the element oriented element 
co-ordinates are selected in the usual way and, starting from 
the harmonic analysis of the angles, the equation of motion 
of each element is derived. The element co-ordinates are also 
considered as system co-ordinates and the equation of motion 
for the whole mechanism is derived in these system co-ordinates. 
The boundary conditions are introduced in the form of a set 
of compatibility equations. USText, the displacements are 
expressed in Fourier series, the constraints expressed by the 
compatibility equations are eliminated and finally the 
differential equation of motion is reduced to a set of 
algebraic equations. Using this set of equations, the 
steady-state deflections and stresses are calculated in 
Section 4.3. Finally, the same set of equations are 
utilised in Section 4.4 to find out the various instability 
zones for the operating speeds. 



4.2 Equation of Motion 


4.2.1 Element analysis 

Eor the method of analysis described in this 

chapter, it is assumed that the mechanism is planar with 

uniformly straight links and has only one input crank rotating 

with constant speed id . The method can be easily extended to 

other mechanisms provided the associated rigid body harmonic 

analysis is already known. 

In case of any element E* shown in Eig. 3.1, 

Fourier series expressions for sine and cosine of the angle 

©g, , are given by eqs . (2.85) and (2.86). Similar expressions 

2 

for the angular speed ug, and the angular acceleration 

ecg, are given by eqs. (2.89), (2.90) and (2.91) respectively. 
Corresponding to the element E' , the element co-ordinate 
axes x-y and the element oriented element co-ordinates H-i 
to ‘tig are selected in the same manner of Chapter III and are 
shown in Eig. 3.1. The element oriented element force vector 
p due to the rigid body inertia forces and other exterml 
forces (external forces having a fundamental frequency 
different from the crank speed will not be considered in 
this chapter) is given by eq. (2.100). 

The series expressions for the rigid body pin 
force Pjj, wg, and obtained from the rigid body analysis 
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(Appendix D), eq. (2.90) and (2,93) respectively, are put in 
the expression of the element oriented element geometric 

Q. 

stiffness matrix It (Appendix B). The resulting series 
expression for k may be written in the following form; 
h. 


■ G 




d 


^ *li 2 

I [£. cos ( 0 - 1 ) on: + k sin Ow-r] 

=1 


h. 


4- io 


I [k cos(o-1 ) 0 ) 1 ? + l:;: sin o o)t] 
0=1 ^ ^ 


(4.1) 


where and are 6x6 matrices. The o)-dependent 

c 

part of E is separately treated for the purpose of the 
stability analysis outlined in Section 4.4. 

Upon substitution of the series expressions for 

p 

0 ^,, o)|, , k , ttg, and p in eq. (3.41), the equation of motion 
for the elem-ent assumes the following form (neglecting p^ 
and n) ; 


mu + 2 U^u + (k° + o)^E^)tI = p 


(4.2) 


= I [ cl cos (o-l) sin 0“''^-] 


0=1 


-1 -a 

n . = VC. m 
0 ^ 


where , 



12S 


-2 


vs. HI 

cl 

h. 


1 

I cos ( 3 - 1 ) (oT + sin j a)T] 

j=1 J J 

h. 


•V 




k5 




D 


J cos( 3 -i) uT+ liT sin 

1=1 


0 


k + ic] 




j — 2j •••} h.i 


(4.3) 


£5 


k'^ 


—3 — •-& 

K< - nc .m + tc .m 

3 3 3 

kt - ns .ffi + ts -ffl®' 

3 3 3 


j 3 — 1j li.] 


and m, k and p are given by eqs. (2.24), (2.23) and ( 2 * 100 ) 

respectively, vc., vs., nc , ns., tc . , ts . are defined in 

J J 0 0 3 J 

Section 2.7 and ^is defined in eq. (3-35). 

Equations similar to eq. (4.2) are derived for all 
elements of the mechanism following the above procedure and 
it completes the element analysis. 


4 , 2.2 Assembly of elements 

Unlike the previous chapters, here the system 
co-ordinates selected are in the directions of the element 
co-ordinates associated with the respective nodes. This 
implies that three system co-ordinates are defined at each 
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intermediate mode within a link in the directions matching 
with the element oriented element co-ordinate system attached 
to that link. Por nodes connecting the end points of severa-l 
links, three system co-ordinates for each of the links are 
defined as done in case of the intermediate nodes and thus 
the total number of the system co-ordinates for such nodes 
is equal to three times the number of the links attached to 
that node. Examples for this kind of the system co-ordinates 
are shown in Fig. 4.1. 

However, this special choice of the system 
co-ordinates does not preclude the use of the code-system of 
assembly described in Section 2.4. On the contrary, the 
effectiveness of the code-system is more apparent in this 
case. After performing the assembly of the element matrices 
m, c, k*^, and the vector p in the manner described in 
Section 2.4, the corresponding system natrices are formed and 
the equation of motion for the whole system is rewritten : 
in n numbers (say) of the system co-ordinates as follows; 

MU + 2C‘^U + 2 C^U + (K° + = 

where, 

o'! = ^ 


P 


(4.4) 




FIG 44 ElEMEMT ORIfiNTED SYSTEM COORDINATES 
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P 


^ [cl cos(3-l)£ax+ sin Iut] 
j=1 ^ ^ 

2, [K. cos(3-1 ) (0-^+ sin j ut] 

j=1 ^ 

V * 3 ■4" 

L [K:^ COS ( j- 1 ) ujT + Kt sin juir] 

j=1 ^ 

hi 

^ i p 

I [P. cos(3-1) a)T4- P. Sin 3b) x] 

3=1 'J *J 


( 4 - 5 ) 


and U = nxl displacement vector in the system co-ordinates 

The nxn system matrices Cl, 0 ^, k!, k 1 , and nxl system 

0 J J J J 3 

12 

force vectors P . , P. are formed from the assembly of the 6x6 

3 3 

-1 -2 —5 -2 -6 --7 

element oriented element matrices c., c., ki, k., k., k. 

3’ 3’ 3’ 3 3’ 3 

12 

and 6x1 element vectors "p., p. respectively (3 = 1, h- ) . 

J J « 

cl 

The matrix C^, taken as proportional to M, represents the 

material damping of the system. The c-oefficient 5 '. may be 

3 

given suitable values for various values of 3'. Prom the 
rigid body force analysis, the input torque maintaining the 
'static’ equilibrium for the rigid body forces is determined 
in series form and is added to P appropriately. The 'speeial 
conditions' of the mechanism are also incorporated in eq, 

(4.4) in the usual way. 
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4.2.3 Constraint equations 

The system co-ordinates selected in the previous 
section do not satisfy the boundaiy conditions existing at 
the joints of the links. The displacement boundary conditions 
at the joints are the compatibility of the displacements. 

Once these displacement compatibility conditions are satisfied, 
the compatibility conditions for the internal forces at these 
joints are automatically satisfied. At the hinged joints, 
the compatibility conditions are; the displacements of each 
link attached to the joint should be same in any two perpe- 
ndicular directions. Since the displacements at the 
supported ends are zero, the corresponding boundary (or 
compatibility) conditions are automatically satisfied if no 
translational system co-ordinates are defined at these ends. 

It is evident that these compatibility equations depend 
only on the angles between the links. 

The set of c numbers (say) of the compatibility 
equations for the whole mechanism may be written in the 
following matrix form: 

SU = 0 (4.6) 



(4.7) 



12 . 

and S^, S^, j = 1, ..♦» are known cxn matrices . 

Explicit expressions for eq. (4.7) are shown in Appendix E 
for two simple mechanisms. 

Equation (4.6) expresses the constraints which 
exist among the system displacement co-ordinates and 
consequently, the equation of motion given by eq. (4.4) is to 
be solved together with eq. (4.6). For this purpose, a set 
of c numbers of the dependent co-ordinates B are chosen from 
the total n numbers of the system co-ordinates. The numbering 
of the system co-ordinates should be done in such a way 
that these dependent co-ordinates form its last part. The 
corresponding partitioned form (the subscript A corresponds 
to the set of the independent co-ordinates) of eq. (4.6) is 
shown in eq. (4.8). 


CSa 


1 


0 

(4.8) 



L°bJ 




or, Uj 

= ^C^A 




(4.9) 

where = - 





(4.10) 


Any set of c nimbers of the co-ordinates from 
the set of the system co-ordinates may be chosen as the 
dependent co-ordinates provided that the sub— matrix 



becomes non-singular (because its inverse is used in eq, 

(4 •TO)). Furtheimore , by writing the constraint equations 
suitably, the sub-matrix can be reduced to a unit matrix. 

The system co-ordinates are now transformed into 
a set of generalised co-ordinates (total number is n-c) in 
accordance with the following co-ordinate transformation; 


U 






(4.11) 


where 



(4.12) 


Using the contragradient law of transformation, the equation 
of motion is next expressed in the generalized co-ordinate 
system as shown by eq. (4.13). 


M*Uq + 2 C*U^v + K*Uq = . P* (4.13) 

where , 

M* = sJmS^ 

C* = S^(c^ + G®)Sp, 

^ ® (4.14) 

K* = S^(E° + 

P* = S^P 

U„ = (n-c )x1 displacement vector in the new general- 

\j 

ised co-ordinate system. 


and 



135 


Since depends on time, all the four co- 
efficients M*, C*, K* and P* in eq. (4.14) are time dependent. 
Therefore, eq. (4.13) belongs to a special class of the 
ordinary differential equations, i.e., the differential 
equations with periodic co-efficients. 

Mathematical theories of the linear differential 
equations with periodic co-efficients rnay be found in the 
standard references^ ” (detailed analysis of a very 
special case of this class of equations, i.e., Mathieu 
equation is available in Refs. (70,(71)). The method of 
the characteristic exponents the averaging method^"^^^ 
and a combination^ of the averaging method with the 

variation of parameters^ ^ have been applied to find the 
steady-state response and the stability of this type of 
differential equations. However, when compared with the 
linear differential equations with periodic co-efficients 
treated in the literature, eq. (4.13) exhibits the following 
special features; 

i) the co-efficient matrix of the acceleration term (i.e., 
the mass matrix) is time dependent, 
ii) the elastic stiffness matrix is also time dependent 
and therefore the natural frequencies of the system 
vary with time. 



iii) the periodic co-efficients of the velocity and 
displacement terms are not small, 

iv) the periodic co-efficients of the velocity term is 

antisymmetric with strong off-diagonal (coupling) terms, 

v) the co-efficient matrices of the displacement and 
velocity term and the force vector contain the 
fundamental frequency of excitation (unknown in case 
of stability analysis). 

vi) the fundamental frequency of the parametric excitation 

is equal to that of the force vector. 

In the context of the above characteristics, it 

appears that the universal method of harmonic balance 

is ideally suited for finding the steady-state response as 

well as the instability zones of eq. (4.13). This method 

has already found extensive applications in the dynamic 

(52) 

analysis of structural systems , and in the study of 
mechanisms with elastic couplings'’ \ The description 
of this method, as utilised in this chapter, is given in the 
following sections, 

4.2.4 Displacements expressed in harmonics 

Since the force-vector and all other co-efficient 
matrices of eq. (4.13) are in the form of Fourier series 
with fundamental frequency « , the steady-state response of 
ecL . (4.1 3 ) will also vary periodically with a fundamental 




frequency u . Consequently, the displacement vector U in 
eq. ( 4.4 ) is an unknown function of time only and may be 
expressed in Fourier series. 



of the harmonics of U and is the total number of sine 
terms taken in the series of eq, (4.15). Tke choice of h2 
depends on the nature of the problem and the accuracy 
desired. In steady-state response analysis, h^ should be 
made equal to the total number of the dominant cosine/sine 
terms present in the force vector P. 

In view of the assumed expression of U, given by 
eq. (4.15), U will be first substituted in eqs. (4.4) and 
(4.6) and the constraints will be eliminated thereafter. 

This procedure is more advantageous than performing the 
matrix operations involved in eq. (4.13) and can be easily 
programmed. 

Upon differentiation of eq. (4.15) with respect 
to time, the expressions for the velocity U and the acceler- 
ation U are obtained below. 

^2 

U = y 3 <0 (- u] sin 3 u'C.+ U? cos juT^) (4.16) 

3=1 ^ 

ff = I - 3^u)^(ul cos 3 toT + u? sin Jojt) 


(4.17) 



Substituting eqs. (4.15), (4-16) and (4-17) for 
if, U and U respectively in eq. (4.4), multiplying the 
series involved and equating the co-efficients of similar 
cosine and sine terms on both sides of the resulting equation, 
the following matrix equation is obtained. 

(- + K^)U® = P® (4.18) 


where , 

S p 

M = general mass matrix with to as co-efficients, 
s 

K = general stiffness matrix whose terms do not 

contain oj. 
s 

P = general force vector. 

TJ® = general displacement vector comprising of the 

12 1 
elements of U. and U., j = 1, h^, and U . 

3 3 ^ o 

(The vectors and matrices in eq. (4.18) are of the order 
( 2 h 2 + 1)n). 

Arrangement of the rows of the matrices and of 

s 

the unknown elements in the \inknown vector U may be done 
depending on the n-umerical method used in the solution of 
eq. (4.18). 

A similar substitution of eq. (4.15) for U in 
eq. (4.6) yields the expanded constraint equation as shown 


below; 
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where S® is a ( 2 h 2 + 1)c x ( 2 h 2 + 1 )ia general constraint 
matrix. The arrangement of the rows and columns of S® must 
follow the similar arrangement of K® and in eq. (4.18) . 

The dependent co-ordinates B (c in numbers), 
chosen during the partitioning of eq. (4.6), is now expanded 
to ( 2 h 2 + 1)c numbers of dependent co-ordinates in the new 
system. Denoting them by B and the remaining (2h2 + l)(n - c) 
numbers of independent co-ordinates by A, eq. (4.19) is re- 
written below in partitioned form; 

isl s|] 

or, U® = S U| (4.21) 

o -'i o 

where, S = - (Sg) (4.22) 

In order to eliminate the constraints (gi-ven by 
eq. ( 4 . 19 )) from the equation of motion (now expressed by 
eq. (4. 18)}, the following transformation of co-ordinates is 
defined: 

U® = (4.23) 

where, 

S* = r-1 (4.24) 

and = ( 2 h 2 ^^^+ - c)x1 displacement vector in the 

■nresent generalised co—oi^inate system. 


{U|} 


iup 


( 4 . 20 ) 


= Co > 
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By the contragradient law of transformation, the mass matrix 
the stiffness matrix K® and the force vector P®, all 
being of order(2h2+1)(n-<3), the new generalised co-ordinate 
system are given below by eqs. (4.25), (4,^6) and (4.27) 
respectively. 


# = 

S*%®S* 

= “aA + “IbS + 

+ kIjS) 

(4.25) 

= 

S*'^K®S* 

= Ka * Kf/ + 


(4.26) 

pS = 

S*'^P® = 

Pj + S^=P| 


(4,27) 


where ^AA' ^AB» sub-matrices 

formed by partitioning of the parent matrices. The subscripts 
A and B denote the independent and the dependent co-ordinates 
respectively. Thus, the equation of motion in the new 
generalised co-ordinate system is given by eq. (4.28). 

(- to^M^ + kS)U^ = pS (4.28) 

This equation will be next used to find out the steady-state 
deflections and the stresses of the mechanism. 

4.3 Displacement and Stress Analysis 

Por a given crank speed oj , the two matrices yf 
and K® in eq. (4.28) can be combined into a single matrix A® 
and eq. (4.28) is reduced to 
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aSd« = P« (4.29) 

where A® = _ (4.30) 

Since the rigid body degrees of freedom are not 
eliminated from the stiffness matrices or or K^, they 
are singular. However, the matrix A® is non-singular because 
of the addition of the non-singular matrix - to K^. The 

unknown displacement vector can be easily found out from 
eq,. (4.30) by inverting the matrix as follows; 

-1 

IjS = (A^) pS (4.31) 

The direct solution (shown by eq. (4.31)) of 
by Gaussian elimination^ ^ ^ is time consuming as the size of 
the matrix A^ is fairly large, i.e., of orfiar {h-c)( 2 h 2 + 1 ) 
and the computations in inverting a matrix increase with the 
cube of its size, if the band-width is not considered. Hor 
this reason, the iterative method of successive over 
relaxation (S .0 .R. ) ^^"* ^ is preferred for solving eq. (4.29). 
In the context of eq. (4.29), this method has the following 
advantages over the direct method of solution, given by eq. 
(4.31). 

i) The number of computations per iteration increases 
with the square of the matrix size. 
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ii) Sine© only 1111© fipsis few iianiioiiics of fh.© fransvepse 
displacemen't and rotational co-ordinates are dominant, 
they can he grouped together to form the most 
dominant suh-matrix of the system. It is well-known 
that a proper rearrangement of the dominant elements 
leads to a large reduction in the number of iterations 
required for solving the equation and the above 
arrangement, though immaterial in the direct method, 
improves the iterative scheme of solution. 

iii) If desired, advantage can be taken of the considerable 
number of zeros present in the matrix A^. This is not 
possible in the direct method of solution (because 
does not have a regular band-width) . 

iv) Of each co-ordinate, the first few harmonics only are 
dominant; consequently greater accuracy is not needed 
in the values of the remaining harmonics. By using a 
stricter convergence criterion for the dominant 
co-ordinates and a less stringent criterion for the 
remaining co-ordinates, the number of iterations can 
be reduced further. This freedom in the choice of 
accuracy is completely absent in the direct method. 

After U® is obtained from the solution of eq. 

(4.29), the vector U is found out from eq. (4.23). As U 

1 . 1 , 2 ' 

is formed by the amplitudes 3 = 1, • ••» 
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the displacement vector U and the acceleration vector U for 
the whole system is known from eqs. (4.15) and (4.1?) 
respectively. In the present method of analysis, the system 
co-ordinates are defined in the directions of the respective 

element oriented element co-ordinate s 3 rstems; so the displa- 

•• 

cements u and accelerations along the element oriented 
element co-ordinates of each element are easily found out by 
simply referring to the integer matrix NS. 

These displacements H and accelerations u for all 
elements of the mechanism are found out in the form of a 
series of the harmonic fimctions of ur and therefore can be 
evaluated for any position of the crank by putting the value 
of the crank angle in these series expressions. 

As the element displacements H and accelerations 
n are known for any position of the crank, the calculation 
of the dynamic stresses within the links follows the same 
procedure as outlined in Section 2.8. The stresses in the 
static restrained structure are calculated during the 
formation of the element load vectors p as mentioned in 
Section 2.8. 

It may be mentioned here that the mass matrices 
and stiffness matrices are shown separated all along only 
to analyse the stability of the mechanism (described in the 
next section) . If the steady-state d isplacements and 
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stresses only are desired, the matrices 1°, in eq. (4.2); 
1°, in eq. (4.4); - A®, K® in eq. (4.18) and - 

K® in eq, (4.28) can be added together (as oi is known) 
thereby reducing the matrix operations to a large extent. 

It will be shown at the end of the next section 
how the present method lends itself to find out the frequency 
response of a flexible mechanism at the cost of little 
additional computations. 

4.4 Stab il ity Analy s is 

Under dynamic loading conditions, a rod may 
undergo three types of instabilities: (i) the instability 
due to the bending resonance conditions, (ii) the parametric 
instability due to the time -dependent axial forces and 
(iii) the instability due to the longitudinal resonance 
conditions. In case of mechanisms, the mass and elastic 
stiffness vary with time as evident from eq, (4.13). 
Consequently, the natural frequencies also change with time 
and the bending resonance conditions, encountered in the 
ordinary elastic systems, degenerate into the parametric 
instability conditions. Ihe natural frequencies of the 
longitudinal vibrations are usually very high; so the 
possibility of the existence of the longitudinal resonance 
conditions is remote. 
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From the theory of parametric Tihrations ^ , it 

is known that the instability zones corresponding to the 

critical (fundamental) frequency u)g_ of the axial forces are 

separated from the stability zones by the periodic solutions 

2 ^ , 4-u 

having periods ~ and Almost in all cases the periodic 

a a 

solutions themselves lie within the instability zones. Two 
solutions of identical periods confine the region of 
instability while two solutions having different periods 
bound the region of : iistability. Though unstable frequencies 
corresponding to all types of subharmonic and ultraharmonic 
solutions exist in the parametric vibrations, the harmonic 
solution and the subharmonic solution of order ^ (as 
mentioned above) only are important from practical point of 
view because all other periodic solutions lie within the 
stability zones corresponding to these two particular 
periodic solutions. 

In view of the above, the stability analysis of 
the mechanism reduces to finding out the critical frequencies 
of the axial forces for which the system will have the two 
types of periodic solutions mentioned above. Since the 
fundamental frequency of the rigid body axial forces is 
identical with the speed w of the crank, the critical 
frequencies of the axial forces mply the critical speeds 
■of the input crank. 



i4(; 


To find onf fhe critical frequencies corresponding 

to the harmonic solution having period it is ^sinned that 

0 )' ' » • 

the homogeneous part of the equation of motion eq. (4.4), 

subject to the constraints eq. (4.6), allows a ’sustained 
(8*7 1 

motion’ given by eq. (4.15). Consequently, the desired 
critical frequencies are given by the solution of the 
homogeneous part of eq. (4.28) as shown in eq. (4.32), 

(- + kS)uS = 0 

The eigenvalue problem of eq. (4.32) may be 
solved by the standard methods^^"^^ . However, in view of the 
singularity of non-symmetry of and the large size of 
these two matrices, an efficient method is described in 
Appendix P. 

The value of hg to be chosen in eq. (4.15) depends 
on the number of instability zones to be determined corres- 
ponding to the harmonic solution of the mechanism. The 
effect of damping, always present in a physical system, 

reduces the width of the instability zones of higher order. 

c 

However, due to the presence of w, the matrices 2wG and 

will contain large off-diagonal terms and the influence 
of the higher harmonic terms on the lower ones may not be 
insignificant. For this reason, the value of h2 should be 
chosen slightly higher than the number of instability zones 
to be determined corresponding to eq. (4.15). 
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To find out the periodic solutions with period 
— , the follov/ing Fourier series expression of U is assumed 
in place of e<l* (4.1 5). 

2h2-1 

U = I (ul cos + U? sin (4.33) 

3 = 1 , 3, 5 3232 

Using this series of U, an equation similar to eq. (4,52) 
is arrived at (following the same procedure outlined before) 
and is given by eq. (4.34). 

(- + K®‘ )U^’ = 0 (4.34) 

The frequencies obtained from the solution of eq. 

4fE 

(4.34) correspond to the periodic solutions with period -jj-. 

The choice of h 2 in eq, (4.33) should be governed by the 
factors mentioned above. 

If the value of h 2 in eq. (4.15) is so chosen that 
it includes the dominant harmonics in the force vector P of 
eq. (4.4), the determination of the critical frequencies of 
the harmonic oscillations may also yield the frequency iresponse 
of the mechanism at the cost of little additional computations. 
If $ is the modal matrix, obtained from the eigenvalue 
solution of eq. (4.32), the following normal co-ordinate 
transformation may be used. 


(4.35) 



I'lS 


wh.6i?6 x\ "th.© disp lac 812161115 tocIjop in the new nopnial 

co-ordinate system* As before ^ use of the above transformation 
in eq* (4.28) yields eq* (4.36)* 

(- a?M^ + K^) n = (4.56) 

where $ and $ (4.37) 

For a particular value of <u , the solution of eq. 
(4.36) is given by eq. (4.37). 

n = (- + K^)"^ (4.38) 

Substitution of this n in eq. (4.35) gives the 
displacement vector in the generalised co-ordinate system. 
The determination of the various harmonics of the elenent 
co-ordinates from the vector follows the same procedure as 
outlined in Section 4.3. 

Once $ is determined from eq. (4.32), the solution 
of eq. (4.36) is rather simple as the order of the matrices 
are small. Since the matrices and z’^ are same for all 
values of the cranZ speed, solution of eq. (4.36) for various 
values of the crank speed requires only the inversion of a 
small order matrix in eq. (4.38) for each value of the crank 
speed (the vector contains the crank speed, in most 

cases, only as a multiplication factor). If the effects of 
the terms due to the Coriolis components of acceleration 
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and the additional tangential components of acceleration are 
neglected, the procedure becomes particularly simple because 
the matrices M and K become diagonal in that case. 

The advantages of the method of analysis described 
in this chapter over the previous method of analysis described 
in Chapter III are; 

(i) It directly gives the steady state solution of a 

flexible mechanism requiring computer time much less 
than the other method 

(ii) It enables to determine the critical speeds of the 
crank for the instabilities due to the subharmonic, 
harmonic and superharmonic oscillations of all orders. 

It is not possible to determine these critical speeds 
from the previous method. 

(iii) If desired, frequency response of a flexible mechanism 
for a wide range of the crank speed may be obtained 
with a computation time less by several orders of 
magnitude than that required by the previous method. 

(iv) In the previous method, the motion during each interval 
of time is described by a few modes of vibration which 
remain valid for that interval only. These modes continu- 
chan-'^o from one interval to another and the 
effects of these varying modes of vibration may be 
pronounced with large external forces or with high 
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( 52 ) 

input speeds ' . While the previous method can not 

take care of the effects of the varying modes, this 
situation does not arise at all in the present method 
as it analyses the mechanism for a general configuration 
(v) In the previous method, the rigid body axial forces 

are considered as time independent during each interval 
of time which destroy the 'frequency modulation' 
effect of the axial forces to some extent. Ihe present 
method treats the rigid body axial forces in their 
actual form. 

(vi) The present method enables to identify the dominant 
harmonics of the displacements which is not possible 
with the previous method. 

The disadvantages of the present method, when 
compared with the previous method, are the following: 

(i) It is applicable to a restricted class of mechanisms 

for which the rigid body harmonic analysis is available, 

(ii) It can not study the transient behaviour of the 

mechanism, particularly at the start of the mechanism, 

(iii) The effects of the elastic axial forces can not be 
studied by the, present method. 

(iv) It may require the pheripheral storage facilities of 

a computer (depending on the order of the matrices ) ^ 
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CHAPTES Y 

EESULTS km DISCUSSIOIS 

5.1 Introduction 

To assess the merits and effectiveness of the 
various analytical treatments presented in Chapters II, III 
and lY, a number of niomerical problems have been solved. 

Eor this purpose, several computer programs have been written 
in a sufficiently general way to accommodate various types 
of mechanisms. In addition, several cases, viz., effect of 
the number of divisions of the links, contributions of the 
various acceleration terms etc. may also be studied by 
controlling +he input data only. These programs were run on 
IBM 7044 and IBM 370/155 computers to solve the numerical 
examples. In some of these numerical examples, the results 
obtained from the present work have been compared with those 
of the previous works. Only limited number of results are 
presented here because of space limitation. 

The following mechanisms are chosen for the 
numerical examples: 

Mechanism 1 (Big. C-1) 

length of linlx 2 (assumed rigid) = _ 

Link 3 (connecting rod): Length ~ J *^ " > cross'^ection = 0.5” 
X 1.0" (larger dimension in the plane of bending), sp. wt. = 
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0.284 Lbs. /in. Young’s modulus = 30 X 10^ Lbs. /in. ^ 

Wt. of the piston = 8.1792 Lbs.| gas force on the piston and 
damping are considered absent. 

Mechanism 2 (fig. 0-1) 

Length of link 2 (assumed rigid) = 1.2"; 

Link 3 (connecting rod): Length = 12", cross-section = 

0.5" X 1.0" (longer dimension in the plane of bending), 
sp. wt. = 0.284 Lbs./in.^, Young's modulus = 30 X 10^ Lbs/in. ^ 
Wt. of the piston = 0,852 Lbs.; damping ratio (C^) = 0.2 for 
the first mode and 0.8 for the higher modes; piston force 
acting on the mechanism is given by eq. (D-4) of Appendix L. 

Mechanism 3 (fig. 0-2) 

Length of the four links are 36", 12", 36" and 30" respectively- 
cross-section of each link = 1.0" X 1.0"; link material is 

■z 

Aluminium with sp. wt. = 0.101 Lbs. /in. and Young’s modulus = 

3 X 10 Lbs. /in. ; moment of inertia of the flywheel at the 

2 

end of link 4 * 0.7 Lb. in. sec ; damping ratio = 0.05 for 
each mode. 

Balancing input torque at the crank end (obtained from eq. 
(2.57)) is assumed to be present at the crank end. 

In the following, the bending stresses calculated 
at the extreme fibre of the links (due to the three types of 

factors mentioned in Section 2.8) are denoted by Oj " 
^^^bending ^^^bendiug respectively. 
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Por the efficient representation of the various 
cases, the followir^ table is used: 


Table 5 . 1 


Case Left side of eq. (3.42) 


A 


II 

= K® 

I 

= K' 

B 

Qa 


II 

= 0 

C 

^a 

II 


i 

= K 

D 

^a 



= 0 

E 



= 0 


E 

T 

K 

= 0 




G All terms are present 


Eight side of eq. (3-42) 
quasi-static 

-Do- 

harmonic series (h^ = 20) 

-Do- 

-Do- 

-Do- 

-Do- 


H All terms are present but -Do- 

the axial forces are 
derived from the elastic 
analysis (described in 
Section 3-4) 


5.2 Basic Analysis 

In this section, results are presented to study 
some of the basic features of the dynamics of flexible 
mechanisms. Eesults have been presented for the first 180 
degrees of the crank rotation in most of the cases. This 
may not yield a complete picture but tte effectiveness 



of the vs.i'ious factors may be assessed fo some extienl; from 
the transient behaviour. 

fig. 5.1 presents a comparison between the stresses 
at the middle point of the follower of mechanism 3 , 
obtained by considering the crank to be instantaneously 
clamped at the support end (it will be hereafter referred to 
as structure) and by following the method of eliminating 
the rr’.gid body degrees of freedom as discussed in Section 
2.5.2. It is seen from the figure that the difference in 
the results is appreciable and the assumption of a clamped 
crank does not yield the correct results, from the figure, 
it is apparent that the frequency of oscillation is reduced 
(with an apparent reduction in stiffness). The reason for 
the reduction in the amplitudes of the vibration may be 
attributed to the fact that the balancing torque applied at 
the crank shaft is taken up by the claaped' end and does not 
contribute to the deformations of the mechanism. The effect 
of such assumption on the steady state solutions is shown in 
Section 5.4. 

fig. 5.2 and fig. 5.3 show the values of 
at the middle point of the follower and the coupler 
respectively of mechanism 5 for various numbers of divisions 
of the follower and the coupler (the symbol i-j-k represents 
i, 3 and k numbers of divisions of the crank, coupler 
follower respectively). It appears from the figures that 
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the stresses (and obviously the elastic displacements also) 
approach their true values when the number of divisions is 
three. Further increase in the number of divisions improves 
the results only marginally though the amount of computations 
increases considerably. 

As explained in Section 2.5.3, an. input torque 
acting on the crankshaft is considered to balance the rigid 
body inertia forces at every configuration. The nature of 
the variation of this input torque with the crank rotation 
is shown in Fig. 5.4 for mechanism 5. The variation of the 
torque at the rocker end is also shown in the same figure. 

It is clear from this figure that the fluctuation in the 
balancing input torque is quite large ani, to obtain a 
reasonably uniform crank speed, a flywheel of considerable 
moment of inertia is needed. The effect of adding a flyv>;heel 
to the crankshaft on the angular rotation of the support end 
of the crank and on at the middle point of the 

follower are shown in Fig. 5.5a and Fig. 5.5b respectively. 
This figures depict the sensitivity of the deformation 
characteristics with respect to the moment of inertia of the 
flywheel (denoted by in the figure) at the crank end. 

Fig. 5,5a shows that, though the angle of rotation becomes 
increasingly smaller with the increase in the values of I^, 
it nevertheless requires a larg© value of to reduce the 
angle of rotation to zero. However, Fig. 5.5b shows that 
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the transition from ’mechanism' to 'structure' is very 

rapidly accomplished by gradually increasing the value of I . 

c 

This is due to the fact that the input torque is largely 
responsible for the vibrations of the crank and its effect 
on the follower is less dominant. 

It is well-known that the natiaral frequencies of 
a mechanism with elastic members depend on its instantaneous 
configurations. Since the configuration of the mechanism 
changes continuously, it is interesting to observe the 
variations of the natural frequencies with the crank rotation. 
As an example, the first natural frequency of mechanism 3 for 
various crank positions is shown in Pig. 5.6. The first 
natural frequency, when the mechanism is reduced to a structure 
by clamping the crank end, is also shown in the same figure. 

It is seen that the first natural frequency is reduced by 
such modification. This fact is also borne out from Pig. 5.1 
as mentioned earlier. It is interesting to note that no 
difference is noticed between the frequencies of the 
mechanism and the corresponding structure when the crank and 
the coupler lies along the same line (i.e., when the crank 
angle is 38 deg. and 238 deg.). It has been further obse^^ 
from the obtained results (not presented here ) that the 
higher order f recjuencies are more drastically reduced when 
the mechanism is conYerted into a 
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It was mentioned in Section 2.8 that the total 
stress (bending, shear or axial) within a link should be 
calculated from three considerations. These three types of 
stresses (for bending only) at the middle point of the 
follower are shown in Tig. 5.7 for mechanism 3 with no 
diTision of the links. It is seen that the stresses 
c^jbending^ calculated from the displacements, are dominant 
over the other two. The observation made in Section 2.8 


that the stresses cTj and cTjjj become small compared to 
when the links are subdivided is substantiated from the 


computed results. At higher speeds, as found in case of 
mechanism 1, the relative importance of Cj and cTjjj increases. 
In case of mechanism 1, with 24000 EPM of the crank, the 
stresses cJj-,- and found to be of the same order 

when no divisions of the links are taken. When computations 


are carried out with three divisions of the connecting rod, 
o^bending i^^^^bending middle point of the connecting 

rod are about average) of at the same 
point. For four divisions of the connecting rod, this ratio 
reduces to about 


5.3 Effects of the Dynamic Factors 


This section is devoted to the discussions of 
the results indicating the effects of the dynamic factors 
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not included in the basic analysis. Again, only the initial 
transient portions of the characteristics are presented to 
study these effects, 

at the middle point of the follower of 
mechanism 5 (with and without the effects of the ri^id body 
axial forces) are shown in Fig. 5.8 and Fig. 5.9 for two 
different speeds. From both the figures, it is seen that the 
effects of the rigid body inertia forces are quite appreciable. 
The rigid body axial forces increase or decrease the 
deflections of the links, depending mainly on the condition 
whether they are in the state of compression or tension. 

The contribution of the rigid body axial forces 
is again shown in Fig. 5.10 for the middle point of the 
coupler at 2400 EPM. It appears that the effect of the rigid 
body axial forces is more pronounced in case of the coupler. 
This may be due to the larger slenderness ratio of the 
coupler which makes it more susceptible to the axial forces. 

In these cases, the rigid body inertia forces acting on the 
instantaneous structure are considered quasi-static (i.e. , 
independent of time during the interval At). If the forces 
are considered to be dynamic in nature in the interval 
(as discussed in Section 2.7), the results are slightly 
modified as indicated in Fig. 5.10. 

Fig. 5.11 to Fig, 5.15 show at the 

middle point of the coupler for various speeds considering 














Fte S.IO COUPLER STRESS FOR CASES A . B % D 

















CA£TQRS _AT' SO-0 



















173 


"tils dynsmic f sctops S6ps.i*Q,'t©ly , Uticugh, iij is s-pps,!*©!!!; fpoin 
the figures that the effects of the additional terms in eq. 
(3.42) due to the centrifugal, Coriolis and tangential 
accelerations increase with the speed, their orders of 
magnitude are quite small. The only other significant contri- 
buting factor is the consideration of the elastic axial 
forces as discussed in Section 3.4. 

When the speed of a mechanism is increased, not 
only the magnitudes of the displacements change, but the 
nature of the variations with the crank rotation also 
undergoes some alterations. 5'ig. 5.14 to fig. 5.17 shov; the 
variations of the stresses for different cases. The apparent 
large difference in the nature of the variations is due to 
the fact that the curves are plotted against crank rotation. 

If plotted against time, the difference would not be very 
large . 

To show the effectiveness of the finite element 
methods considering all the aspects discussed in Chapters II 
and III, the problem of a slidercrank mechanism (mechanism -1) 
has been solved. The deflection at the middle point of the 
connecting rod is shown in Fig. 5* 18. Comparing this ■with 
the similar results presented in References 4 and 5, it 
is seen that the nature of the variation of the deflections 
at the middle point of the connecting rod is in good 

..... '■ f 

Q.gTBBTiion'fc with those ohtaiiied xb the two works mentioned 
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The discrepancies may be attributed to the difference in the 
methods of analysis. Transverse displacement of a point (at 
a distance of 4.8" from the crank end) of the connecting rod 
of mechanism 1 is again plotted in Fig. 5.I9 at the same 
speed, but with case A. The importance of the different 
dynamic factors may be appreciated by comparing Fig. 5. 18 
with Fig. 5.19. It appears that when the dynamic factors 
are ignored, the nature of the vibratory motion is completely 
changed. This is due to the fact that the mechanism is 
running at a very high speed and the dynamic factors, 
particularly the rigid body axial forces, become more pronounced. 

5.4 Direct Steady State Solution 

Since the finite element technique is very 
versatile and powerful, a new approach for the direct steady 
state solution using the finite element techniques has been 
proposed in Chapter 11 . In this section, the numerical 

results are presented and compared with those from the j 

■■ ' ' ' ' ' ^ ^ " ■■■■■ ■ ■' ' , 

previous works in some cases. Both the slider crank and | 

four— bar mechanisms have been solved using this approach. 

Solving mechanism 3 With this method , the i 

rotation Ug (in radians) at the coupler-end of the follower 
is found in the f ollowirig form: 




r CriftH CTI-' C^THE COMNECI roo 
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Ug = 10 ^(0.0253 + 8.12 00802 + 3.31 cos2©2 + 2.29 C 0 S 3 Q 2 

+ 0.0481 cos 4©2 - 1.59 cos6©2 +9.86 sin©2 + 7.50 sin202 
+ 5.81 sin3©2 sin4©2 sin5©2 sin602) 

where ©2 is the crank: angle. 

Similar expressions are also obtained for the displacenEnts 
at other nodal points of the links. Fig. 5.20 to Fig. 5.22 
show the variation of at the middle of the crank, 

follower and the coupler. Sadler and Sandor^"'*^^ have solved 
the same problem following the lumped parameter approach. 

The nature of the variation of the stresses will not be same 
as that of the deflections because the rigid body displacements 
do not take part in strainirig the member. However, the 
difference is found to be less pronounced in case of the 
follower. The nature of variation of the stress at the 
middle point of the follower shown in Fig. 5.22 is almost 
identical with that of the deflection at the same point found 
by Sadler and Sandor^”**^^ . The same problem has been solved 
considering the crank to be instantaneo\;ely clamped and the 
results are shown in Fig. 5.21 and Fig. 5.22. 

Mechanism 1 has also been solved by this approach 
and the deflection of the middle point of the connecting rod 
for one complete cycle is shown in Fig. 5.23. Comparing 
this With Fig. 5 . 18 , becomes clear that the steady state 
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pa'fc'fcGm can b© coinple'tely differenl; from fhe transient: part. 

If the steady state solution is determined following the 
conventional methods (with two divisions of the connecting 
rod), almost 15 minutes of computer time is needed in IBM 
7044- in contrast to 2 minutes taken by the new method of 
analysis. In complicated mechanisms where the number of 
co-ordinates is high, the computation time is reduced by an 
order of magnitude if this new method is adopted instead of 
the conventional finite-element approach. 

Mechanism 2 has been solved for various values of 
the crank speeds and the noimalised deflections of the 
middle point of the connecting rod (for one complete cycle) 
are shown in Fig. 5.24. The normalisation has been done to 
reduce the difference in the order of magnitude of the 
deflections for various crank speeds. Since tte rigid body 
inertia forces increase with the square of the crank speed, 
the normalisation has been done by dividing the actual 
deflection by where = angular velocity of the crank/ 
first circular natural frequency of the connecting rod, 
considering it to be a simply supported beam. Though ; 

mechanism 2 is not, strictly speaking, a simple degree of | 
freedom system, its characteristics exhibit a good resemblance 
to those of a single degree of freedom ^stem. When the 
amplitude of the deflections at the middle point of the 
connecting rod are plotted against the frequency ratio 
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the curves obtained (Fig. 5.25) resembles the frequency 
rGspotisG cuz'VGS of 9. dsiiiipsd. singlo dGgrss of firGsdom sysl/Gui 
with Gxcitation baing proportioml to tho squaro of tho 
forcing fraquoncy. Pig. 5.24 also indicatas a changa in the 
phase batwaan. the excitation and deflection whan the crank 
speed is varied. Pig. 5.26 shows the variation of phase 
angle with the frequency ratio Prom this figure also it 

is clear that the phase-shift characteristic is not much 

different from that of a damped single degree of freedom system. 

( 5 ) 

Viscomi and Ayre ^ solved the same mechanism for the steady 
state deflections with = 0.3 and found the steady-state 
amplitude to be about 0.72 in. Pollowing the method of 
Chapter IV, the steady-state amplitude is found to be, 0.79 in. 
The study of the speed-response characteristics of mechanisms 
is very much facilitated when the present method is adopted. 
Otherwise, considerable amount of computation time is 
required in the conventional methods. 

In this chapter (except Pig. 5.5) comparisons 
between a 'mechanism' and 'structure' have been made by 
assuming that no flywheel exists at the crank-shaft of the 
mechanism. This has been done only to demonstrate the 
difference between the results obtained from the two 
extreme cases. It is evident that the results for a 
mechanism with a fl 3 mheel at the crank-shaft will lie 
between these two cases as - shown in Pig. 5.5. 
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Only a few numerical examples have been solved in 
this chapter as the major objective of the present work is 
to introduce the met nods of analyses and demonstrate their 
efficacies rather than to study a few specific mechanisms 
exhaustively. 

Since the numerical values of the examples solved 
in the previous works are inP.P.S. units, the same system 
of units is chosen in the examples cited in this chapter to 
facilitate direct comparison. 
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CHAPTER VI 
GONCIUSIONS 

Several important conclusions may be drawn 
following the numerical results presented in Section 2.9 and 
Chapter V. These may be listed as follows; 

(i) When a mechanism is converted into a structure by 
imposing an aritifical constraint (viz., clamping 
the support-end of the crank), the vibration 
characteristics including the natural frequencies 
are altered. The natural frequencies of a mechanism 
are reduced as a consequence of such transformation. 
Therefore it should not be included in the analysis 
of flexible mechanisms. 

(ii) To pbtain even moderately accurate results, it is 
necessary to divide the links into more than one 
element.. The optimum choice of the number of 
divisions of a link depends on its relative flexibility. 

(iii) The stresses obtained from the analysis of ’static’ 
and ’dynamic’ restrained structures are less than 
those obtained from the displacement analysis alone. 

The relative magnitudes of the former with respect 
to the later type of stresses depend on the speed 

of the mechanism. With increase in the number of 
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(iv) 


(v) 


(vi) 




(Yii) 


(viii) 


divisions of a link, the stresses obtained from the 
restrained structures gradually decrease. 

The general method for the rigid body analysis 
described in Section 2.9 takes computer time less 
by an order of magnitude than that taken by the 
existing matrix procedures. 

The effects of the rigid body axial forces on the 
transverse vibrations of the links are appreciable 
and in the kineto-elastodynamic analysis of high- 
speed mechanisms these effects should be taken into 
account. 

The dynamic nature of the rigid body inertia forces 
(instead of considering them to be quasi-static) 
should be taken into consideration for better 
accuracy of the analysis. 

When the axial forces are obtained from elastic 
deformations of the links, their effects on the 
.transverse vibrations of the links are even mo3?e 
pronounced. 

The effects of Coriolis components of accelerations, 
the additional normal and tangential components of 
accelerations on the vibrations of the links are 
not significant, at least in the transient part of 


a f our»bar mechanism. 



Computer time taken by the new method of analysis 
presented in Chapter IV to obtain the steady state 
solution is less by an order of magnitude when 
compared with the conventional methods. 

The total number of harmonies to be taken in the 
new method of analysis depends on the number of 
dominant harmonics present in the force vector. 
S-uf f iciently accurate results have been obtained 
with eight harmonics in a case of a slider crank 
mechanism described in Chapter V. 
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Appendix - A 

Description of Computer Programs 

In order to solve the numerical examples in 
Chapter V, several computer programs have heen written. Of 
them, only two major programs are described below. These 
programs are sufficiently general in the sense that these 
can analyse (i) a mechanism or a structure (having no rigid 
body degrees of freedom), (ii) with arbitrary number of 
divisions of any link, (iii) with or without axial deforma- 
tions of the links, (iv) with or without the effects of 
rotary inertia and shear deflection ,(v) a quasi-static or a 
complete dynamic problem of a mechanism etc. 

Program 1 

The purpose of this program is to solve eq. (3.42) 
for the types of mechanisms cited in Chapter Y. All cases 
described in Chapter II (with the exception of Sections 2.6 
and 2.9) and Chapter III are included in this program. The 
major steps of computations are: 

(i) Read the input data (consisting of geometrical, 
material properties etc. special conditions like 
damping, external loading etc. and the necessary 
codes) . 



(ii) Find, "bh.© miglosj sng'u.lm’ vslociilssj sq^usi*© of "th© 
angular vslocities, angular accelerations and the 
absolute longitudinal and transverse accelerations 
of the left end of each link, all in series of 
harmonics of the crank angle, 

(iii) Find the element oriented element matrices m, £ 
and for each link. 

(iv) Starting with a crank angle, find the quantities 

calculated in step (ii) by summing up all harmonics 
for the current crank angle. 

(v) From the rigid body analysis, find the total rigid 
body axial pin forces acti.ng at the left end of each 
link. For the analysis described in Section 3.4, 
this step may be omitted. 

(vi) Transform the element matrices of eq. (3.41), the 

initial displacement vector and velocity vector into 
the system oriented element co-ordinates by use of 
the rotation matrix R in eq. (2.26) . 

(vii) Find the stresses for all elements in the ’static' 
restrained structure. Also find the element load 
vector and the element geometric stiffness matrix 
for all elements in the system oriented element 


co-ordinates. 



(viii) Assemble the element matrices and Teeters to form 
the system matrices M, C^, K+Ak'^ and and the 
vector P (expressed in harmonics) of eq. (3.42). 

Also form the initial displacement and velocity 
vectors for the whole system from the assembly of 
the corresponding terms found in step (vi). 

(ix) Add ’special conditions' of the mechanism to the 
relevant terms of the system matrices. 

(x) Find the transformation matrix T defined in eg. (2.41). 
Also find the input torque from eq. (2.57) and add 
it to P. 

(xi) Solve the eigenvalue problem of eq. (3.45) to find 
the spectral and the modal matrices for the whole 
system. 

(xii) Transform the system matrices, vectors and ’initial 
conditions' into the normal co-ordinate system. 

(xiii) Solve the equation of motion, eq. (3.42) (using the 
method of Section 3.3). 

(xiv) Transfoim back the displacements, velocities and 

accelerations from no 2 mal co-ordinate system ’obtained 
in step (xiii) to the system co-ordinates . 

(xv) Transforming the above quantities in element oriented 
element co-ordinates, find the stresses and 
within the elements and the nodal displaceiiBnt and 



velocity vectors in element oriented element 
co-ordinates to be used in step (vi)as 'initial 
conditions' for the subsequent crank position. 

(xvi) If the analysis of Section 3.4 is used, find the 
increment in the axial forces of the elements and 
test the convergence criterion. If the convergence 
is achieved, go to the next step. Otherwise, form 
the incremental geometric stiffness matrix, add it 
to the previous stiffness matrix and go back to 
step (xi). 

(xvii) Find the new crank position by adding the prescribed 
increment to the current crank angle and go back 
to step (iv). 

Program II 

The purpose of this program is to solve eq. (4.29) 
for the types of mechanisms cited in Chapter "V. The program 
may be divided into the following major steps; 

(i) Head the input data (consisting of geometrical, 

material properties etc. , 'special conditions' like 

external loading, damping etc, and the necessary 
codes). 

(ii) Find the angles, absolute angular velocities, square 
of the absolute angular velocities, absolute angular 



accolGra-tiions j ijli© absolute longitudinal and, trans- 
verse accelerations of the left end of each link, 
all in series of harmonics of the crank angle. 

(iii) I'ind the element matrices m, c®', k^ and 

(excluding the contributions of the element geonetric 
stiffness matrix) of eq. (4.2) for each link. 

(iv) I'rom the rigid body analysis, find the harmonics of 
the torque at the crank end and the rigid body axial 
pin forces acting at the left end of each link. 

(v) Eind the stresses for all elements in the ’static’ 
restrained structure. Also find the element load 
vector and the element geometric stiffness matrix 
for all elements in the element co-ordinates. 

(vi) Assemble the element matrices and vectors to form 
the system matrices M, and (K^ + u^k"'^) and 
vector P of eq. (4.4). 

(vii) Add 'special conditions' of the mechanism to the 
relevant terms of the system matrices. 

(viii) Find the constraint matrix S of eq. (4.6). 

(ix) Construct the matrix (- wV + K®) and the vector 

P® of eq. (4.18) from eqs. (4.4) and (4. 15 ). to (4.17). 1 

(x) Construct the roatrix S® of eq. ( 4 . 19 ) from eqs. (4,6) 
and (4.15). 

(xi) Eliminate the constraints to form the matrix A^ of | 



(xii) Solve eq. (4.29) for the displacement vector #, 

(xiii) Rearrange tho elements of to find the displacement 
at each co-ordinate, expressed in series of harmonics, 
(xiv) Ror any position of the crank, find the displacements 
(and accelerations) at the element co-ordinates. 

(xv) Rind the stresses cTjj and within the elements. 

The special subroutines used by the above two 
programs are described here very briefly. 


Routine SCSKBY 


In tho above programs I and II , multiplications 

of tho harmonic series of the following kind are required 

very frequently, 
h^ ^2 

I [a, cos(i-1)0 + b. sin i©] I [c. cos(o-1)e + d. sin j©] 

hj 

= I [e-, cos (1-1)9 + f. sin 19] (A-l) 

1=1 

where the scalars h^ , h 2 , b^ (i = 1, h^ ) and c^, 

d . (3 = 1 , ...,h 2 ) are known. Ror a given value of h^, the 
scalars f^ (1 = 1 , . . • , h^) are to be found out from the 

multiplication of the two harmonic series at the left side 
of eq. (A-1 ) . Ror this purpose, a routine- SCSKIW has been 

written which returnsthe coefficients e^, f^ (1 - •••> 

h^) as output after performing the indicated series multi- 
plications (for any value of 9). 



Routine SCALUN 




In step (iv) of Program II, it is necessary to 
solve eq. (A-1 ) for the unknown scalars a., b. fi - i 

1 ^ i ^ ‘ 5 * • • > 

h-, ) when the scalars (3 = 1 , . . . , ) , e^, f^ (1 = 

1, h^ ) and h^ (= hg = h^) are known. A routine SCALOT 

has been written which returns the coefficients e ,, b . 

1 A 

( i = 1 , ... 5 h^ ) as output (for any value of 0) after 
equating the coefficients of the similar harmonic terms on 
both sides of eq. (A-1). 

Routine UMRQY/ !; 

■ -i 
i 

This routine has been written to obtain the ‘ 

matrices M® and K® respectively of eq. -(4.18) as output v/hen ■ 
the matrices M, C^, C°, K°, k"''" of eq. (4.4) and the scalar 
h2 of eq. (4.15) are supplied as input. The same routine 
also returns the 'combined value of (- + K ) of eq._ (4.18 

if the value of u) is supplied. Por the assumed series of U ; 
in eq. (4.15), the output matrices are non-symmetric . i 

However, by using a multiplication factor i suitably, these 
matrices are rendered a symmetric form so that they may be 
used more efficiently in the numerical procedures. Heedless 
to say that the final solution should be modified accordingly 
to obtain the correct results. The special features of the 


routine are; 
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(i) The routine is equally applicable if the unknown 
displacement vector U is expressed in the form of 
eq. (4.33) instead of eq. (4.15). With little 
alterations, the routine may be used for the stability 
analysis of the subharmonic and superharmonic 
oscill<ations of all orders. 

(ii) Thu columns of the output matrices may be separated, 
if dosired, into the following six different groups. 
The olemonts of those groups are associated with 

(a) the preferred co-ordinates (these co-ordinates 
should form the first part of the system 
co-ordinates) and the preferred harmonics (i.e., 
a few pre-selected dominant harmonics) 

(b) the preferred co-ordinates and non-preferred 
harmonics 

(c) the non-preferred co-ordinates (excluding the 
preferred co-ordinates and the dependent 
co-ordinates) and preferred harmonics 

(d) the non-preferred co-ordinates and the non- 
preferred harmonics 

(e) the dependent co-ordinates and the preferred 
harmonics 

(f) the dependent co-ordinates and the non-preferred 


harmonics . 
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The second feature of the routine enables to use 
certain numerical procedures (for example, Gauss-Seidel 
iteration procedure ^ ) very efficiently. As this routine 
uses the peripheral storage facilities of a computer, it 
can handle veiy large matrices. However, in spite of large 
size matrices, the routine takes very little computer time 
as no multiplication is involved (in the real sense). 

Routine UNTSTR 

This routine has been written to find the matrix 
S® of eq. ( 4 , 19 ) from eq. ( 4 . 6 ) by using the assumed harmonic 
expression eq^ (4.15). With the same routine, similar 
matrix maj be obtained corresponding to the assumed relation 
eq. (4.33) for the system displacement U or for other assumed 
relations of U corresponding to the subharmonic and super- 
harmonic oscillations of various orders. This routine is 
very similar to the routine UNEHOW and possesses the same 
two features, like UNKNOW, it also takes very little 
computer time (as no multiplication is involved) and uses 
the pheripheral storage facilities. 



Appendix - B 

Explicit Expressions of the Element Matrices 

The explicit expressions for the element stiffness 
matrix k, element mass matrix m, element force vector p and 
the element geometric stiffness matrix k , given by eqs. 
(2.23)} (2,24)} (2.25) and (3.36) respectively, in the 
element oriented element co-ordinate system are shown in 
eqs. (B-1 ) to (B-4) respectively. 

^Mx (B-4) is the axial force (+ve along +ve 

x-direction) acting at the left-end of the element E' and 
(Dg, is its angular velocity. It may be noted that when the 
links are subdivided into several elements, the matrix 
and the vector p are different for these elements even when 

matrix m^ given by eq. 

antisymmetric 

0 

-(4+3T)l 0 

0 42+40f 0 

0 -(6+5?)l 0 0 

(B- 


the links are uniform. 

The explicit form of the 
(3.35), is shown in eq. (B-5). 


m 


P A1 

T20{^+(f) 


0 


72+40<p 


0 


(6+59)1 0 

0 -(1 8+209) 

I8+2O9 0 

-(4+59)1 0 



1 2EI Symmetric 

6EI1 (4+9)EI1 
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Append 0^ - 

Direct ()n„lcul.-i,t;ic)n of the Transformation Matrix 
Introduction 

A general procedure for calculating the rigid body 
transformation matrix of a mechanism has been outlined in 
Section 2.5.2. For some mechanisms, the transformation 
matrix may also be determined directly from the geometi’y. 

Two illustrative examples are given below for a(plane.) 
slider-cranl-c and a cranlc~rocker mechanism. 

Slider-orank mechanism 

For the v-th division of the link§(cranlc) of the ' 
slider-crank mechanism shown in Fig. C-1 , the co-ordinates 
right-end A in the system co-ordinates X-Y are; 

^ ““sSg ; ^ f®-'' > 

Where, 

Lr, = length of the link 2 

dg = total number of divisions of the link 2 
V = division number of the element 
and Sg = angle made by the link 2 with X-axis. 

Let the rotational co-ordinate of the support-end 
of the crank be chosen as dependent co-ordinate, Correspo.-.dii 
to a small Increment SGg “ angle O^, the yariatione 
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FI& C-1 SLIDER CRAMK MECHAMISM 



and oYjv in the oo-or,Jinat« X^, \ are obtained from eq. 

( 0 - 1 ) as follows; 


6X 


■A 


_ \> 3 x 13162 ^^2 * 


cosOg 6©2 (C-2) 


For the v-th division of the link 3 (connecting 

, . n- r 1 the co-ordinates X^j of its right- 

rod) shown in I'lg. ^he co 0 Bn 

end B in the system oo-ordinates X-Y are. 


L.^ 

cogG^ + ’ 

= Lg sin@2 + 


(C-3) 


where 


I. 


A - total number of divisions of the link 

3 


length of the link 3 


^ -1 4 .V. innir '5 with the X-axis, 

and ©V = angle made by the link ^ 

. .+ 60 in the angle ©2 A 

Corresponding to a small increme 2 ^ 

AY and 6 © in the co-ordinates X^, Yg 

variations fiX., SY^ and 

^ fG-3) as follows; 

the angle ©^ are obtained from 


6X. 


■B 


and 


6Y. 


■B 


• r\ aq — V sin©’z 

sin ©2 662 ^ ^ 

La 

Q AO + ■'^ A ^^3 

L 2 cos ©2 <^5 ^ 


(C-4) 

(C-5) 


I'or the slider-c 


• m the following relation holds . 

X'EdIc inechB-D-isiHj 


sinQ-x 


^ siB.©n 

. «»*• -r / 


(C- 6 ) 



2 ID 


n V 


6©.. 


Iig COS02 


60. 


_ _ — =. ■■■■ ■ . . " ' A— OWo 

3 cos0^ 2 

EeplaoemBnl, oi: 8 O 3 from oqs. (0-4) and (0-5) by eq. (C-7) 

yields ; 


(C-7) 


6X. 


3 


Ip sin0, 

{- 1.2 slnOg + V ^ oosSg) 602 


61- 


3 


L. 

(1, 2 00002 - v ^ oosOg) 602 


(C-8) 

(C-9) 


Thorol-oro, the elemerts of the transformation 
matrix T (dofitiocl in eq. (2.41)) corresponding to the 00 
ordinates Ov^, ”'’3 ^ obtained from 

gq . (C-2 ) as : 


T V . 


V - 1 — sinOp ; Iv 

Up ^ 


=v?cos02; Tv3=1 

2 “ ^ dg ^ ^ 


(C-10) 


Similar elements of the matraoc 1 corresponding ^ 
the oo-ordinates Uv^. Ovg and TI., at the point B are obtaine 


;rom eqs. (C-8), (0-9) and (0-7) as follows: 
Tv. 

I ~ 

I. 


Tj sin©3 
^ ^ COS©o 

Ijg sinOq + cos©^ ^ 


T' 


Tv. 


1-2 ^os©2 - cosSq 

Ii„ GOS0p 
*" 1*3 COS 03 


(0-11) 
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Wiuolvuiigiu 

li'or the crank-rocker mechanism shown in Fig. C-2, 

t. v V of" tVip ri^'ht—eni C of the v — th division 

the co-ordinates Xq, Yq oi the ri^n^ 

of the link 4 (rocker) are; 




^0 = h'-' ' 


(C- 12 ) 


(C-13) 


where ; 


... icnigth of the link 4 


d ^ total number of divisions of the link 4. 
and «! - -'.n-Mo by tho libk 4 with the X-axie. 

,0. a =:n,nh-hoo.er mechanise, the following 

. (40 ) 

relations liold ’ 


;os9. 



(1-0. 

'3-"“""^^ ’ (cp^-2c^ 00302+1 . (C-15) 


3 in©. - ^ 


ucl =oa0,)«ine5 


Cos©/ 


cosQ^ (c ^ T 2 


(C- 16 ) 


A Vi'To ' + Go COS 02 )^ (^ 

oosS,) - Obse, *1-(0i+02 2 


sin©^ = 


sin©^ ( 0^+02 




FIG C-a CRANK- ROCKER MECHANISM 
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Lg sinOp + Lj sinOj + sin©^ = 0 

Lp cosOp + cosOj + cos©^ = 


(C-18) 

(C-T9) 


where 


- Lp - 

2L^ 


»D j Ij Q 

p I ^ ^ — £ 

V O T T 

2 


Ci *ss 




^2 t2 t2 t2 

. 4 * h - ^4 + h 

°4 “ — ~ye^^ 


c < _ ~ 

°5 " L, 


(C-20) 


C 4 = ^ 
“6 


Corresponding to ^ small increment ^©p in ©g, variations 6©^ 
and 6©^ in the angles ©^ and ©^ respectively are obtained 
from eqs. (G--18) and (0^19) as follows: 


60 . 


L- sin(©^ “ 

__2 „ A — , — 6 ©„ 

L, sin(©, - ©4 j 2 


(C-21 ) 


60 


4 


Lo sin(©o - ©•z) 

2 , 4 > 5 © 

sin(©j - © 4 ) 2 


(C- 22 ) 




Substitution of 60.^, obtained from eq. (C-21), in 


eqs. ( 0 - 4 ) ' 

and (0-3) yields 



6Xb = 

(- L2 sin02 - V sin© .7 

3 ^ 

Y' ) 6©2 

(C-23) 

1! 

L, 

(Lg cos©2 + V ^ cos©^ Y' 

) 6©2 

(C-24} 

where Y ' 

Lg sin(6p - 0^) 

"■ *“ I, .5 sin(©' - ©4) 


(C-25) 


Corresponding:^ to a smail increment 60^ in the crank angle ©g, 
the variations SIq and 6Yq in Xq and Yq respectively are 
obtained from eqs. (G~12), (C-13) and (C-22) as follows; 


v 6 X, 


6Y. 


I^(- 1 + f-) sinS^ Y" 602 

4 


L^(- 1 + h) oosS^ Y" 602 

'4' 


(C-26) 


(C-2Y) 


(C-28) 


Ig sin (©2 - © 2 ^ 
where 'lY" = sin('^ 

Thus, the elements of the transformation matrix T 
corresponding to the co-ordinates at the points A, B, and C 
are given by eqs. (C-10), (C-29) ^'nd (C-30) respectively. 


Tv. 


lo sin©o ~v Sin©^ Y' 


Tv2 = '^2 °°®®3 


(C-29) 


Tv, 


Y 



22i 
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Tv2 


L^(-. 1 + GOS©^ Y" 

(G-30) 

T 


Y^t 



The 

arraw^emont of these elements 

in the matrix T 


in aocordanoo w.i.th tl'io numbering of the system co-ordinates 
may bo convoniuntly done by using the integer matrix KS . 



Appendix - D 
Rigid Body Analysis 

Procedures for direct kinematic analysis as well as 
static-force analysis of a slider-crank mechanism and a crank 
rocker inochaninin ax'o outlined below. Since the mechanisms 
are considered to be composed of rigid links only, it is 
convenient to use minimum number of co-ordinates for each 
link, i.e., to treat each hinge-joined link as a truss member 


S 1 id 0 r-c rank Mechanism 

For the slider-crank mechanism shown in Pig. C-1 , 
the absolute angular velocity and angular acceleration 
of the connecting rod are determined by differentiating eq. 


(C-6) with rcispoct to time: 


0.^ = 

cosOn 

2 2 ^ 

],i.y cosO-^ 2 

P 0 


©., .-s 

5 

L,, oinOp 4* sin©^ 

©2 

CO'SO-'r 



(B-D 

(D-2) 


where io the uniform angular velocity of the cr 

. ■ ( 88 ) 

The piston acceleration is 

i OOS 202 + . . 


(a 


, 572 ^ 


Sin ©2^ 




where n 



2?0 


The following analytical expression^5) 
force "E'q (acting on the position) may be assumed 

52 MjLg 


for the fi-as 


M^L5(3.75 - 1.25 cosO^ _ 


2-5 cosO^) (d_ 4 ) 

where is the mass of the link 3 (connect tv. 

2 rod). 

The standard sh£ipe function for a 

^ truss member E' 

(Pig. D-la) is 

S C 5 


a 


- c -I 


0 1-C 0 ^ 

where 5 and 5 are eis defined in eq. (2.I3). 

Using this sliape function in eq. ( 2 . 25 ) the 
statically u’qui.'vaionl modal forces due to the rip-id body 
inertia fore us acting 611 the member E' ares 


(D-5) 


P'1 


' “Sx "1 

P2 




6 

P3 


^Mx ^'^Ux 



%y 


(D-6) 


where lVi|,, , is tliu macs of the member E' and 
and arc as defined in Section 2,3- 


Por the uniformly rotating crank 

0 ana = 


tfe “ -^-My “ '*^Ny 


* 2^2 



9 '■.) 7 

< 4 



pm D-l FORCES, m h SUDER CRANK MECHANISE 
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For till' contu-ct i u,'; roci », 


II 

1 

^0*1^0 C0S(©p 

“ ©3) 5 

Kr 

Nx 

= a. cos©.. 

V ■ “ 

Q^lip win, (©2 

©3) ; 

Kr 

Ry 

(D-8) 

= - sin©^ 

Using eqs 

. (:ih- 7 ) and 

(.D~8) in eq. 

(D-6), 

the nodal forces F.| 

to (shown in Fig. 

B-l ) are obtained. 

From the free-body 

diagrams 

oC bho crank 

(IMg. D-lb) 

and the 

connecting rod 

togc tho.r 

Willi fciuj piutnn (Fig. O-lc), the 

following equations 

result; 

4 :i'' 

3 4 - R^ + R^ 

= 0 



F^ + 1'' 

4 + ”2 + h 

= 0 

(D- 9 ) 


T, .. ( 

R 4 *'4)I'2 

= 0 



I'Y + n 

6 = 0 




and coo©^ = -- Pg. + M^oc^ - (F^ + F^) cos0^ + Fg sin0^ 


where .iu the onbput torque at the crank end and R^ 
are pin roantiono on the links as shown in Fig. B-r1b and c. 
From the corajja t’j 1) i 1 i by ui' the internal forces at the hinged 
joint B (Fig. the following relations are 

" ^ ^ ^ ^ ; (D-10) 

iq. oLn(«g - ©3) - Rg cos (©2 ~ ©3) 

Aftor oJ.irnination of the unknowns Rg, R3, R^ and 
Rg from thu uqu. (h-Q) and (D-10), the desired axial rigid 



• f R , torque are finally obtained 

body pitt iorco.J i 

nr'Ank-rocker.ig acIiailj^^ 

The absolute angular velocities ©3 and ©^ of the 
coupler and the rocker respectively of the mechanism shown 
in Pig. ^-2 obtained in eqs. (D-11 ) and (D-12), after 
dividing eqru (0^21) (0-22) by 6t. 


© 


'3 




'4 




Y"»o 


(D-11 ) 
(D-12) 


* 

, 5ri non (0-25), (0-28) and ©p is the 

where y*» Y ‘ ^ ^ ^ 

uniform angular veiocity oi the crank. 

Upcni differentiation of eqs. (D-11) and (D-12) 

With ro.pect to ti»o, tho abociut. aagalar accelerations §, 
and 9^ of tho coupler and the rooker respectively are gxven 

by eqs. (D-13) nnd (D-14). 




L,,e? cos(6„-0j t ccs(9g- 9^)^hi... 

*■ L3 sin(©3"“©^) 


(D-13) 


008(9^-63) - 1.6? (D-14) 


T‘4 sinCSj-h 


Tho nodal forces to Pi 2 

btained fron eq. (D-6) when the following relations are 


used . 




r>'2 im a rocker, mechahisja 
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For the crank: =s 

For the coupler; 


V = o -na V = -e2i 


2-^2 


(I)-15) 


Mx 


©gLg cosCg^-q^^^ ^ 
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•My 


©^Lg Din (©2 --9^). ^ 


= V ■■ 


For the rocker: A^^ 


♦ O 

T 

^ ^ Xi 


■% ^ S[j + 

V = - 84 h= V = V = 


(D-16) 


0 


From the free~bod^r ^ 

<aiagrains of thp n 

4 {fig. D-Ka to o), the foli„„,. ^ and 

equations are obtained: 

+ itj i 


(D-17) 


- 0 

= 0 

I ^“4 = 0 


«a t t I'g t 


I'l t (H^ + F^)! 

Ej + + Fj + ^ Q 


“b + ^8 

’t ^6 


(D-18) 


0 

0 


'P . 
0 


(^10 ^1G^\ = 0 


where R- to li. are the pin reactirtnc. 4.u -, 

T H ^ '"GXions on the links and 

To are the torquoh at the crank and rocker ends as shown in 

Fig. D-Z. Thu ooinpatibillty equations for the Internal forces 

at the hlnKod jolnt.s may be written from Pig. D-2d and Pig. s-z. 




o 


R oou(« 2-*'^;5) - % 

✓ 

uin(^^.,,-«-5) - \ con(Og-tt-^) 
“ R<|q 

Eg = sin(©3-«4) - H^o 


(D-19) 


Asnuming tlmt is known, the desired pin forces 
R R R and the toi’quu muy be obtained by sol'ving eqs. 

I y sa# 

(J3-1B) and (U-1h). 



233 


Appendix ~ E 

BxpreBu;ionn Tor Harmonic Analysis and Constraint Equations 
giider>»crank mechanism 

Tlsinfr the same symbols of Appendix C, the relevant 

harmonic express Iona, correct up to -2L- a.re: 

n^® ’ 

C0S&2 = ( 1-a|~'5H2-3ba^-.126ag) + (a^+4a2+1 5aj+56a^+210a^) 

C01J2W., (a,,f6a.^+28a4+l20a5)cos iOg + (a^+Sa^+iSa^ ) 

1C)^g)oos802 + oobIOO^ (E— 1 ) 

“ ©I 12Ba^+512a^) + (2a^--16a^-102.4a^- 

312ai^)cDa20^ - (Ba2+32aj+102.4a^+292.5714a^)cos4©2 





Constraint ediiationsi. 


Prom U.n comna U,M 11 fcy of the rt in placements at the 


hinged Joint H h-ld) ntui at the piston powitton C (Pig 

D-1e)^ tile i'oJ 1 uwi.iig »'»huitlnnw i'»)nult,s 
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in eq. 


sin©2 

COS©.z 


expressions for the fractional terms 

(B~ 11) are: 

- n ^ ^ •5a^ + 15a2+87.5a^+401a^)sin©2 - (0.5a^+7.5a2 

+ 52.5a^+294a^)sin3e2 + d .5a2+17.5a5+126a^)sin5©2 


(2. 5a^+31 .5a^^)sin792 + 3.5a. sin99p ] 


(E-12) 


cosO, 


1+o..jH-9a2+50a^+245a^+1134a^) - (a.j + 1 2a2+75a^+392a^- 
1890a^)cos2G2 + (3a2+30a3+196a^+1080a )cos4e - 


L3+56a^+405a3)cos602 + (7a, + i90a )cos8Q, 


9^5 coslOOp 


(B-15) 


where 


a-j *fco are as defined in eq. (E-5). 


Crank -roc leer mechanism 

The following series expressions result from 
the theory of complex variables; 


COS 92 ^3 


I (-04)^“^ cos iGp), 
i=1 


(c.p 


+ 20^ cos©2 + 1 


s inO^ 


(Gf < 1) 


(C^ ) + 20^ cos©2 + 1 


I sin i©p 

i=1 ^ ^ 


(E-14) 


* ^ 
3 ^ 





vj ,i!KJC -^-rTTT* and ■ 


14(5 
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less than unity, the 


follfwi n,*’; r;f|nntl(jn i;; obtained from binomial expansions. 


Vi . (CM 4 a;, coo(i,,)‘ 

I 'Am 


whore 


= (1+ ^ COSO 2 ) (1+ qrr^osQ^) 

= Vl-cj^ l”* d]_d^ cos^^3-2 


(E-I5) 


d . 


sr. 1 


dl * 

JU. 




ri. 

2 1 


^2 

up to i~1 ternB) ( grlj ) 


i > 1 
(E-1 5 ^ 


i-1 

up to j*1 terms) (q i_y)‘'^ > J "* 


and h-j in an dufincd in eq. (2.87)* 

lJ5jin#c the following trigonometrical identity 


(89) 


■iH-l „ ...n,, _ QQQ n6 + n cos(n-2)© + ^ cos(n-4)© + ... 


eoe « 


is cos© 


for n = odd 


|(n~l) i(n+1) 


(E-1 7) 




r. 


for n = even 


B) 

2' 


in eq. (E-lb), the following Fourier series expression 


IS 


obtained: 



2:;S 


a-l 

= .1 cos(j-i)02 

3 

3 . h. are known functions of C’ and C^. 

where df , 3 ='»•*• M 

Bse of e^s. (B-H) and (B-18) in e^. (O-U) to 
(C.17) and .altiplloations of tne series yield tne harmonic 

expressions(^°hor the angles In the following forsas^ 

Yi 

f [d"': cos(3-l)©2 


(E-1B) 


COSO-Z — 

^ 3=1 


h. 


sin© 7, 


0 


cos© , 


"1 

i 

3=1 


I [d® 000(3-1)02 t '^3 2^ 


(E-19) 


h, 


l" [d^ 003 ( 3 - 1)62 ■*■ 


3=1 

h 


3,n0, = .< [d- 003(3-1)02 ^ 

t 3=1 

i 11.1 h. are known oo -efficients, 

where d^ to d. , 3 = 1. • ' ' 'll 

Oonstsih^^auatio^^ 

displacements at the hanged 

■ , ! .p.2d and l-2e) yieldsthe followxng 

joints B and C (Bio* 

equations*. 


Ih 


(a Q ) - Bn sin(©2"®3^ *” ^ 

U^_ cos(© 2“®3^ dg 2 



0 


2.19 


sin(Q 2 -^e„) _ u, cos(©^-0,) = 

4 5 ^ ^6 “^3 


U 


j - U, cos (9^-0.) - U, 

^7 M 4' '"d.. sin(0,-- 


(E-20) 


= 0 


U- 


^dg + sinCQ^^©^) - cos(©^-0^) = o 

where to are nodal displacements as shown in Pig. 

* I; 0 


D-2. Considering j S'Hd as displacements at 

3 4 7 8 

the dependent co-ordinates, the constaint equation for a crank- 
rocker mechanism, obtained from eq, (E-20), is written in 
eq. (E-21) (shown just after eq. (E-11)). 

The harmonic series expressions for the rigid body 
axial forces and the balancing torque at the crank end may be 
obtained from the rigid body analysis of these two mechanisms 
as described in Appendix D, For this purpose, the special 
routine SCSKMW, described in Appendix A, has to be extensively 
used for the multiplications of the harmonic series involved 
in the calculations. In the process of eliminating tte 
unknowns, it is further necessary to use another special 
routine SCALUN described in Appendix A (for example, to find 
out the harmonic series expression for the reaction from 
the last equation of eq, (D-9)). 

It may be noted that it is also necessary to use 


the same routine SCSKHW repeatedly in the analysis procedures 
described in Section 2. 7» - in Subsections 4.2.1 , eqs. 

(B-1 1 ) and (E-21 ) . 




Appendix - P 

Solution of large Eigenvalue Problems 


An efficient iterative procedure for large 
generalised eigenvalue problems has been outlined in Ref. 
(90) which uses generalised Rita method along with subspace 
iterations. Without going into the theories involved, the 
necessary steps of the procedure are given below. 

■ let the undamped generalised eigenvalue problem be 
given by eq. (E-1 ) 


K 0 = M$ 


(P-1) 


where, 

K = stiffness matrix of size n X n (say) 

M = mass matrix of size n X n 

q2 = m X m a diagonal spectral matrix whose principal 
diagonals are sqtiares of the lowest eigen- 

frequencies (m numbers) 

and <5 = n X m modal matrix whose columns are the eigen- 

vectors of the system. 

Starting with an initial modal matrix the 
iteration scheme is; 


i) Find by solving eq. (P-2) 


K$. = 


(P-2) 


(i = 1 , 2 , 5 , » • • • ) 
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j_i) Calculate 


*^1 


i 


M. = 

1 i 1 




ill) Solve lor the following reanoed elgensystem 


KlQ^ = Vi“i 


(M) 


(where Q. 1= the modal matrix corresponding to the 


^ — 

reduced eigenvalue problem in eq. (J-4))to find 


, 2 
and 


iv) Find an improved approximation to the eig-nv 


$iQi 


(M) 


-ic. -hchieveds otherwise replace 
v) Stop if convergence is achi 

i l,y i+1 and go bach to step (l). 

,or emcient calculation of step (i), the matrix 

matrices by Oholeskl's metho 
K is decomposed in two triangular matnce 

■tion^®''! Moreover, the eigenvalue problo 
of decomposition eigenvalues and 

(V ll is to he solved for all eig 

defined by eq. 

eigenvectors. For this purpos , 

. . d jaoobi' B method . Bitishauser- s method etc . ) 

g,enorallBed J mo+Vind 

, to the iterative Power method. 

suitable compd-r a good modal matrix 

of choosing a go 

«.e procedure 

.nd a suitable convergence oritorxo 
and therefore is not repeated here. 


It may be noted that, by virtue of the application 
of Ritz's analysis the iterations converge very rapidly 

(usually within eight iterations) for the lowest freQuoncies 
if the eigenvectors are orthogonal. When K and M are both 
symmetric, the eigenvectors are mutually orthogonal with 
respect to K or M. However, following the method of analysis 
presented in Chapter IV, the mass matrix # in eq, (4.32) 
becomes asymmetric and the stiffness matrix K® in the same 
equation, though still remains symmetric, is singular due to 
the presence of the rigid body degrees of freedom. In view 
of these special features of the matrices and the 
above iterative procedure should be implemented along with 
the following modifications. 

Assuming a single degree of freedom mechanise, the 
total number of the dependent rows (or columns) in of eq. 
( 4 . 32 ) is w^ = 2h2 + 1 , h 2 being already defined in eq. (4.15). 
Let eq. (4.32) be partitioned in the following way. 


“aa 

1 


^A 

+ 

^AA 



cil* 

— 1 

M§a 

L BA 



uS 

i 

^ba 

^BB 




Denoting the number of elements present in U ox 

eq. (4.32) by n^ = ( 2 h 2 + 1 ) (n-c ) number of elements 

in uf and u| are equal to n 2 = n^ - and w^ respectively 
(hg, n and c are defined in Chapter IV). The matrices 



mS k® etc. are the partioned blocks of and # whose 
dimensions ensure the matrix maltlplioations in eq. (J- • 

sollowlng the procedure described in Section 2.5.2, the 
eigenvalue problem, defined in eq. (4.52), may be reduced to 

eq. (^'-T). 

(, eV + r')Tr' = 0 (J-T) 


where , 


= Kf 




X 


* 

g 


'AA 

m|a “Ib' 

T; = - 
= of - 

a.- vnci+TviY 1 is defined by 
and the rigid body transformation matn Jg 


(?- 8 ) 


the following relation; 


Kfj,Tg + K|B = 0 


(P-9) 


vrX - gq (P-7) is non=-singnlar and 

The matrix K m eq, ^ 

>,e decomposed into two triangular matrrce . 
thcrei-orc may , iterative procedure perfoxms 

^■•-WK-c' the above mentioned xtera P eiaen- 

^ „ -nf eigenvalues and eigen 

■Pnr the m numbers vOf 

iterations lor . rewritten in the form 

vootors slitiulihdeously, eq. 


of eq. (!'-■') follow = 


^ 


fi 
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where ^ is the 112 X m modal matrix and is the m X m 
(diagonal) spectral matrix as defined in eq. (P-l ) . 

Since the matrix 1/F in eq. . (I'-10) is asymmetric, 
the eigenvectors in the modal matrix $ may not he orthogonal 
with respect to or Kp", Consequently, the convergence may 
he slow if the iterative procedure is directly applied to 
eq. (P-IO). To overcome this difficulty, the iteration 
procedure should he carried out in two different phases. 

In the first phase, the iteration may he started 
hy considering only the upper half of the matrix I/P" and then 
making its lower half equal to this upper half. This modi- 
fication renders the matrix iP' a symmetric form. Next, 
steps (i) to (v) of the iteration procedure should he carried 
out until a moderate convergence criterion is satisfied. 
Since the factors contributing to the asymmetry of are 
small compared to the basic ssrmmetric property of #, it in 
expected that the eigenvalues and eigenvectors thus obtained 
are fairly close to the desired actual quantities. 

The second phase of iterations begins with the 
actual form of the matrix 1? in eq. (I’-TO) together with the 
eigenvalues and eigenvectors obtained at the end of the 
first phase of iterations. Instead of using suhspace 
iterations, the second phase consists of inverse iteration 
in the original space of dimension n 2 and evaluation of 



Rayleigh’ s quotients in each iteration. The relevant steps 
of the second phase of iterations are; 

(vi) Find by solving eq. (R-ll ) 

vii) Calculate Rayleigh’s quotients 


"i.j 


. K 0. . 

lil 

. M 0. . 

IjD i»D 


{¥-^2) 


viii) 


p 

where $- - and nf . denote the j-th column of and 

i > i] ^ » D 

j-th diagonal of fi? respectively, j = 1, m. 

Stop if convergence is achieved; otherwise replace i 
by i+1 and go back to step (vi). 

ix) Find the modal matrix for eq. (F-6) by applying the 
last relation of eq. (F-8) to the modal matrix 
obtained at the end of the second phase of iteration-, 
(an equation similar to eq. (2.53) may be easily 

deduced for this purpose). 

While the rapidity of the convergence rate in the 
j-lfUt ptase of iteration largely depends on the orthogonality 
ol- the modes obtained from step (Hi), similar rapidity in 
(,v,o second phase depends only on the good estimates of 
and (it does not regulre any orthogonality between the 
modes for rapid convergence). As a good estimate of and 



R. is obtained from the first phase of iteration, the inverse 

iteration in step (vi) converges very quickly (usually in 3 to 

4 iterations ). The faster rate of convergence in step (vi) 

is further enhanced by the quadratic rate of convergence of 

the method of Rayleigh's quotients in step (vii). 

The size of the matrices and if being fairly 

large, periphoial storage facilities of a computer nay have 

to bo used if the high speed storage capacity is very limited. 

In tho.t case, the change-over of I?" from symmetry to asymmetry 

in the two phases of iteration does not demand any extra 

location in the core memory. . 

In the context of the specialities of the matrices 

involved, the above scheme appears to be more efficient than 

the popular practice adopted for asymmetric matrices, i.e., 

successive applications of Power method, inverse iteration 

(91 ) 

of Wielandt and deflation procedure^ . 



